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ABSTRACT 


The  BOXPLX  of  a user  specified  cost  function  is  evaluated 
as  a control  system  design  method.  The  emphasis  is  on  compen- 
sator design  for  multiple-input-multiple-output  plants,  but  the 
technique  should  be  applicable  to  any  control  system  with 
adjustable  parameters. 

The  engineer  must  define  the  plant  dynamics,  inputs  and 
outputs,  the  desired  output  responses  for  the  specified  inputs, 
weighting  factors  for  the  performance  measure,  the  type  of  com- 
pensation to  be  used,  the  free  parameters,  and  initial  values 
and  constraints  on  the  free  parameters . 

The  program  gives  a direct  search  of  the  optimum  values, 
within  the  specified  constraints,  for  the  free  parameters. 
Evaluation  of  the  results  is  given  in  graphical  plots  of  the 
output  time  responses  including  the  desired  output  response  as 
well  as  the  compensated  response. 

The  Transfer  J-Iatrix  Method  is  used  to  provide  an  analytical 
check  on  the  accuracy  of  the  method  and  the  procedure  is  illus- 
trated with  a two- input -two-output  system  example. 
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I . INTRODUCTION 


As  control  systems  become  more  and  more  complex,  design 
engineers  are  faced  with  an  increasingly  difficult  problem  in 
the  development  of  compensators.  Single-input-single-output 
linear  systems  can  be  compensated  using  Bode  plots  or  root 
locus  plots  and  trial  and  error  methods  to  meet  the  system 
specifications  [1] . When  compensating  multiple-input-multiple- 
output  linear  syscems,  use  can  be  made  of  transfer  matrix 
techniques  described  by  Ogata  [2]  to  achieve  total  decoupling 
and  the  desired  output  responses.  Since  decoupling  during 
the  transient  period  of  real  systems  may  require  controls 
which  exceed  physical  constraints  imposed  by  the  system,  it  is 
sometimes  desirable  to  decouple  during  steady  state  conditions 
and  allow  coupling  during  the  transien-t  period  [3]  . 

Non-linear,  single-input,  single-ouput  systems,  wherein  the 
non-linearity  can  be  lumped  into  one  element,  can  be  analyzed 
and  compensated  effectively  using  the  phase  plane  method  (for 
second  order  systems)  and  describing  function  analysis  [2,4]  . 

Some  of  these  Ji^thods  require  extensive  mathematical  manipulations 
and  become  rather  unwieldy,  but  do  produce  the  desired  results. 

In  all  of  the  analysis  and  design  approaches  mentioned  above, 
any  one  of  several  well  established  simulation  programs  can  be 
used  to  evaluate  the  design. 

To  reduce  the  tedium  involved  in  the  design  of  compensators, 
frequency  domain  and  time  domain  computer  optimization  schemes 
have  been  developed,  which  will  set  free  parameters  in  the 
compensator  to  yield  a system  response  which  most  nearly 
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approximates  a reference  response  specified  by  the  user.  Lima 
used  this  approach  to  address  the  problem  of  controller  design 
for  ships  engaged  in  underway  replenishment  tsl.  MacNamara 
applied  this  method  to  an  autopilot  design  [61,  and  Vines 
developed  a generalized  program  for  compensator  optimization 
in  single- input-single-output  systems  [7]. 

The  following  chapters  present  a generalized  discussion  of 
transfer  matrix  analysis  and  compensation  of  multivariable 
systems,  optimization  techniques,  and  performance  measures. 

A two-input-two-output  system  is  then  compensated  using  transfer 
matrix  techniques,  and  using  parameter  optimization.  The 

/ 

results  are  then  compared. 
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II.  MULTIVARIABLE  SYSTEM  COMPENSATION 

A.  TRANSFER  J-IATRIX- ANALYTICAL  METHODS 
1.  Introduction 

A linear,  multiple-input-multiple-output  or  multivariable 
system  can  be  described  by  a transfer  matrix,  which  is  simply 
an  extension  of  the  transfer  function  concept.  Consider  the 
two- input- two-output  open  loop  plant: 


Figure  2-1  Multivariable  Plant 


The  Block  Diagram  of  the  plant  can  be  represented  in  the 
general  sense  as  shown  in  Figure  2-2. 


Figure  2-2  Generalized  .Multivariable  Plant 
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Algebraic  equations  for  the  outputs  in  the  independent  variable 
S,  can  be  written  directly: 


C1(S)  = Gpll(s)Rl(s)  + Gpl2{s)  R2(s) 
C2(s)  = Gp21(s)Rl(s)  + Gp22(s)  R2(s) 


If  the  inputs  Rj (s)  and  the  outputs  C^(s)  are  each  considered 
as  vectors  of  one  dimension  (i,  j = 1,  2)  the  equations  can  be 


written  in  matrix  form: 


Cl(s)” 

_C2(s)^ 


GplKs) 
Gp21 (s) 


Gpl2{s) 
Gp22 (s) 


Rl(s) 

R2(s) 


The  number  of  inputs  and  outputs  can  be  increased,  and  need  not 
be  equal  so  that  in  the  general  case: 


CKs) 

C2(s) 


C^(s) 


GplKs)  Gpl2(s) Gpij(s)f  [r1(s 


Gp21 (s) 


Gpil (s) 


Gpij (s) 


R2(s) 


C(s)  = Gp(s)R(s) 


where  Gp(s)  is  the  transfer  matrix  (ixj)  of  the  plant. 


2 . Feedback  Compensation 

Transfer  matrices  can  be  used  to  describe  the  compen- 
sation of  multivariable  systems  as  well  as  the  plant.  Consider 
the  two- input- two-ouput  plant  above  with  feedback  compensation 
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as  shown  in  Figure  2-3. 


Figure  2-3 

Multivariable  Plant  with  Feedback  Compensation 


Algebraic 

equations 

for  R^(s)  and 

R^Cs)  can  be  written: 

Rl(s) 

= UKs) 

H ll(s)  CKs) 

- H12(s) 

C2(s) 

(6) 

R2  (s) 

= U2(s)  - 

H 21(s)  CKs) 

- H22(s) 

C2(s) 

(7) 

with  C,  R,  and  U two  element,  one  dimension  vectors  these 
equations  can  be  rewritten: 
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or  in  the  general  case: 

R(s)  = U(s)  - H{s)C(s)  (9) 

where  H(s)  is  the  transfer  matrix  (jxi)  of  the  compensation. 

From  the  uncompensated  plant: 

C(s)  = G (s)R(s)  (10) 

Substituting  equation  (9)  for  R(s)  in  Equation  (10) : 

C(s)  = Gp(s)  [ U(s)  - H(s)  C(s)]  (11) 

Rearranging : 

-[I  + G (s)H(s)]C(s)  = G (s)U(s)  (12) 

Since  Gp(s)  is  an  ixj  matrix,  H{s)  is  a jxi  matrix  the  product, 
Gp(s)H(s),  is  a square,  ixi  matrix.  Assuming  that  [I  + Gp(s)H(s)] 
is  non-singular  the  inverse  can  be  taken. 

Premultiplying  both  sides  of  equation  (12)  by  [I+G  (s)H(s)]”^: 

^ P "• 

C(3)  = [I  + G (s)H(s)]"^G  (s)U(s)  (13) 

C(s)  = G(s)U(s)  (14) 

The  compensated  system  transfer  matrix  (ixj)  is  then: 

G(s)  = [I  + G^(s)H(s) ] '^G^(s)  (15) 

- ~ , p ~ ..  p 

The  plant  with  feedback  compensation  can  be  depicted  more 
clearly  as  shown  in  Fig'ure  2-4. 
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Figxire  2-4 

Matrix/Vector  Representation  of  Multivariable 
Plant  with  Feedback  Compensation 


To  determine  the  feedback  compensation,  H(s),  required  to 
achieve  design  specifications,  the  dynamics  of  the  plant  must 
be  known  and  put  into  transfer  matrix  form,  G^(s).  The  desired 
closed  loop  dynamics  must  then  be  reflected  in  the  closed  loop 
transfer  matrix,  G(s).  Then  through  matrix  manipulation  the 
compensation,  H(s),  can  be  computed  in  closed  form  provided 
that  G(s)  and  G^Cs)  are  square  and  non-singular. 

Premultiplying  both  sides  of  equation  (15)  by  [l+G  (s)H{s)] 

~ -P 

and  then  postmultiplying  by  G(s)~^: 


I + G (s)H(s)  = G (s)G(s) 


-1 


Rearranging : 


G (s)H(s)  = G (s)G(s) 

^ -•  -•  y -N. 


-1 


- I 


Then  multiplying  both  sides  of  equation  (16)  by  G (s)~^ 

P 

H(s)  = G(s)"^  - 5 (s)“^ 


(16) 


(17) 


(18; 


3 . Cascade  Compensation 

Cascade  compensation  can  be  treated  in  a similar  fashion. 
Consider  the  same  basic  plant  with  cascade  compensation  as  shown 
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Multivariable  Plant  with  Cascade  Compensation 


The  algebraic  equations  for  R^(s)  and  R2(s)  are: 

Rl(s)  » Gcil(3)U^(s)  + Gcl2(s)U2(s) 

R2(s)  = Gc21(s)Ul(s)  + Gc22(s)U2(s) 


"rICsH  Tg 

Jl2(s)j  [G 


Gcll(s) 


Gc21{s) 


Gcl2 (s) 
Gc22  (s) 


1 pJKs)! 

J L^2(s)J 


in  the  general  case: 


G(3) 


Gp ( 3 ) Gc ( S ) U ( 3 ) 
Gp ( 3 ) Gc { 3 ) 


If  Gp(3)  i3  a 3quare  non-3ingular  matrix  the  compensation 


Gc(3)  is  given  by: 


Gc(s)  = Gp 


The  compensated  plant  can  be  depicted  as  shown  in  Figure  2-6. 


Figure  2-6 

Matrix/Vector  Representation  of  Multivariable 
Plant  with  Cascade  Compensation 


4.  Cascade  and  Feedback  Compensation. 

Cascade  and  feedbac.k  compensation  can  also  be  combined 
as  shown  in  Figure  2-7. 


Figure  2-7 

Matrix/Vector  Representation  of  Multivariable 
Plant  with  Cascade  and  Feedback  Compensation 
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k 


E(s)  - R(s)  - H{s)C(3)  (25) 

0(s)  = G (s)E(s)  (26) 

C(s)  = Gp(s)U(s)  = Gp(3)Gc(s)  [R(3)  -H(3)C(3)]  (27) 

C(s)  » Gp(s)  Gc  (3)  R(3)  - Gp(3)Gc(s)H(s)C(3)  (28) 

Gp(3)Gc  (3)  R{s)  » [I  + Gp(s)Gc(3)H(s)  ]C(3)  (29) 

If  Gp(3)  is  an  ixj  matrix  Gc(s)  will  be  a j x j matrix,  and 
H(s)  will  be  jxi,  in  which  case  Gp (s) Gc (s) H (s)  is  a square 
inatrix  (ixi).  If  [I  + Gp(s)  Gc(s)H(s)  ] is  non-singular,  it  can  be 


inverted  and: 

C(s)  =»  [I  + Gp(3)Gc(s)H(s)  ]"^3p(s)Gc(s)R  (s)  (30) 

C(3)  =*  G(s)  R(s) 

The  compensated  system  transfer  matrix  G(s)  (ixj)  is: 

G(s)  * [I  + Gp(3)GC(s)H(s)  ]"^  Gp(s)Gc(s)  (31) 

Equation  (31)  contains  two  unknown  matrices,  H(s)  and  Gc(s). 

One  of  these  must  be  selected  arbitrarily  and  then  in  certain 
special  cases  the  other  can  be  solved  explicitly.  For  example: 

“ [*  0 

In  which  case: 

G(s)  - [I  + Gp(s)Gc(3)  Gp(s)Gc(s)  (33) 

[I  + Gp(3)Gc(s)  ]G(s)  - Gp(s)Gc(s)  (34) 

G(s)  = Gp(s)Gc(s)  [I  - G(s)]  (35) 


If  Gp(3)  is  a square  matrix  G(3)  and  Gc(s)  will  be  square.  If 
Gp(s)  [G(3)  ^ - I],  and  G(s)  are  non-singular: 
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Gp (s) 

= Gc(s)  [G(s)"^  - I] 

(36) 

Gc  (s) 

= Gp(s)“^  [5(s)"^  - I]"^ 

(37) 

If  Gc(s)  and  {[G(s)”^  - I]Gp(s)}  are  assiimed  to  be  non-singular 
equation  (37)  can  be  rearranged  for  convenience; 


Gp ( s ) Gc ( s ) 

= [G(s)"^  - I]"^ 

(38) 

Gp  ( s ) * 

[G(s)"^  - I]"^  Gc(s)"^ 

(39) 

Gc(s)"^  = 

[G(s)~^  - I]  Gp(s) 

(40) 

Gc  (s)  = 

{ [G(s)"^  - I]  Gp(s) 

(41) 

This  result  can  also  be  obtained  directly  from  equation  (37) 
using  the  matrix  identity  (AB)”^  = A~^. 

In  each  of  the  systems  discussed  above  the  compensation 
cannot  be  solved  for  in  closed  form,  except  in  special  cases 
as  indicated. 
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B.  PARAMETER  OPTIMIZATION  - COMPUTER  METHOD  | 

1.  Optimization  Technique 

The  Complex  Method  of  constrained  optimization  (BoxPLX)  ^ 

was  utilized  in  this  program.  The  algorithm  suggested  by  | 

Box  (9)  was  modified  by  the  programmer,  R.  R.  Hilleary,  Naval 
Postgraduate  School,  to  find  the  constrained  minimum,  and  j 

includes  an  integer  programming  option,  as  suggested  in  (10)  . i 

The  Complex  Method  is  a modification  of  the  unconstrained,  direct 

I 

search.  Simplex  Method  introduced  in  (11) . Direct  search  methods 
compare  function  and  constraint  values  only  in  searching  for  ^ 

i 

the  minimum  value  of  the  function.  These  direct  search  methods  j 

have  been  widely  used  because  of  their  robustness,  reliability, 
and  ease  in  programming  and  use.  They  have  widespread  applic-  j 

ability.  The  one  drawback  is  that  they  are  generally  less  j 

efficient  than  gradiant-based  techniques  (12)  . A more  (.'etailed  ! 

I 

discussion  of  the  Complex  Method  is  included  in  the  program 
documentation.  I 

In  this  program,  the  main  program  simulates  the  reference  J 

I 

output  and  then  calls  FUNCTION  BOXPLX  which  in  turn  calls  FUNCTION  j 

FE.  FUNCTION  FE  calls  SUBROUTINE  Plant  which  simulates  the 
compensated  system.  FE  then  computes  the  performance  measure  I 

I 

and  returns  to  BOXPLX.  BOXPLX  keeps  track  of  the  compensator  j 

1 

parameter  values,  computes  a new  set  of  free  parameters  and 
calls  FE  again.  This  iteration  continues  until  termination  as 
described  in  the  documentation.  See  Appendix  A for  more  ! 

information  on  the  program. 
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2.  Performance  Measures 

In  order  for  a numerical  method  for  parauneter  opti- 
mization to  produce  a result,  a decision  process  must  be  com- 
pleted, the  result  of  which  is  the  set  of  parameters  which  is 
defined  as  "optimum."  This  term  optimum  is  an  enigima  of  sorts. 
The  program  addressed  here  utilizes  a comparison  of  time  domain 
response  of  a compensated  system  with  a reference  response 
specified  by  the  user.  The  word  optimvm  above  is  dependent  on 
the  designer's  knowledge  of  what  reference  response  is  best 
suited  to  the  output  he  is  trying  to  control.  Furthermore, 
with  the  exception  of  very  simple  linear  systems,  it  is  diffi- 
cult to  match  a reference  response  exactly,  so  the  designer 
must  decide  if  the  comparison  at  steady  state  is  more  or  less 
important  than  that  during  the  transient  period.  Additionally 
very  high  frequencies  with  periods  less  than  one  half  the 
integration  step  size  may  be  filtered  out  by  the  digital  process 
in  the  computer,  and  not  be  present  in  the  simulated  output. 

The  designer  is  faced  with  some  decisions  which  he  must  make 
objectively  before  the  "optimum"  set  of  parameters  can  be 
computed. 

The  decision  process  mentioned  above  is  based  upon  a 
performance  measure  which  is  an  indicator  of  the  "goodness" 
of  a given  set  of  parameters  being  varied  in  the  optimization 
process.  The  designer  must  select  the  performance  measure;  the 
performance  measure  is  then  minimized  by  numerical  methods, 
since  the  objective  here  is  to  have  the  system  response  match 
the  reference  response. 
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The  term  selectivity  is  used  to  describe  a quality  of  the 
performance  measures.  If  three  dimensional  space  is  considered 
with  the  performance  measure  plotted  as  a function  of  two  free 
variables  in  a compensator  of  the  system,  the  resultant  topology 
should  include  one  or  more  minima  for  the  performance  measure. 
The  minima  with  the  lowest  value  are  referred  to  as  the  global 
minima,  although  the  minimization  performed  is  constrained. 
Selectivity  is  directly  related  to  the  steepness  of  the  slope 
of  the  contour  in  the  area  immediately  adjacent  to  the  minima. 
This  selectivity  is  a function  of  the  performance  measure  used 
as  well  as  the  plant  and  compensator  dynamics  (2,8)  . 

Several  performance  measures  have  been  suggested  in  numerous 
sources,  some  of  which  will  be  discussed  here.  In  each  case 
the  performance  measure  will  be  represented  by  the  variable  J. 

Consider  a single- input  single-output  system  with  output, 
c(t),  and  reference  response,  r(t).  The  error,  e(t),  for  this 
case  is  defined  as  the  difference  between  the  reference  res- 
ponse and  the  system  response  or: 

e(t)  = r (t)  - c{t)  (42) 


Integral  square-error.  The  defining  equation  for  finite 
time  is : 


J 


e^ (t) dt 


(43) 


This  performance  measure  is  not  very  selective  and  tends  to 
emphasize  large  errors  and  ignore  small  ones  (2,8). 
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This  performance  measure  is  more  selective  than  those  mentioned 
previously  and,  as  in  the  case  with  the  integral  of  time  multi- 
plied square-error  criterion,  this  performance  measure  tends  to 
emphasize  errors  late  in  the  transient  period  and  in  steady 
state  and  deemphasize  those  which  occur  early  in  the  transient 
period  (2,8) . 
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The  optimization  of  the  performance  measure  can  be  addressed 
in  two  categories;  the  case  where  the  reference  response  speci- 
fied is  the  forcing  function;  or  the  reference  response  is  some 
function,  usually  a second  order  response,  which  the  designer 
believes  to  be  a model  response  for  his  application.  In  the 
first  caise  the  peak  overshoot,  setting  time,  rise  time,  and 
natural  frequency,  are  determined  by  the  dynamics  of  the  com- 
pensated system  and  by  the  performance  measxire  selected.  How- 
ever, if  the  designer  desires  a second  order  response  with  a 
certain  minimal  rise  time,  he  can  specify  a second  order 
reference  response  which  meets  this  requirement,  and  use  the 
integral  absolute  error  or  integral  error-squared  performance 
criteria  to  optimize  the  compensation  parameters. 

Provisions  for  both  of  these  approaches  have  been  included 
in  the  program  in  that  the  user  can  reference  a second  order 
response,  where  he  specifies  the  damping  factor,  natural  frequency 
and  gain,  or  he  can  utilize  the  table  loop-up  feature  in  which 
ha  can  specify 'any  reference  response  he  desires. 

Other  performance  measures  can  be  developed  using  classical 
measures  of  system  performance  such  as  maximum  overshoot,  time 
for  the  error  to  reach  its  first  zero,  time  to  reach  maximum 
overshoot,  settling  time,  frequency  of  oscillation  of  the 
transient  or  steady  state  error.  However,  performance  measures 
based  on  these  characteristics  would  be  complex  and  increase 
computation  time,  and  their  desirability  over  those  previously 
discussed  is  questionable. 

The  performance  measure  used  in  the  optimization  program 


(see  Appendix)  is  the  sura  of  the  weighted  performance  measures 
of  each  output  of  the  multivariable  system. 


M 

J = E (W^)  (47) 

i=l 

The  siabscript  i denotes  the  output  for  which  J-  ^ is  calculated. 

1.  j K 

The  subscript  k denotes  one  of  four  performance  measures  which 
can  be  selected  by  the  user  to  be  applied  to  each  output  in 
turn : 

Integral  square  error 

Ji,l  = Z (Rj,i  - C.  .)^  (48) 

j=i 

Integral  of  time  multiplied  square  error 

N 

Ji,2  = Z ((Rj,i  - C.  J^)(iAT)  (49) 

j=l 

Integral  of  absolute-error 

N 

Ji,3  * E I Rj,i  - C . I (50) 

j=l 

Integral  of  time  multiplied  square-error 

N 

Ji,4  = Z I Rj ,i  - C . J (iAT)  (51) 

j ' 

j = l 

N represents  the  number  of  integration  steps  in  the  simulation, 

AT  represents  the  integration  step  size,  M is  the  number  of 
outputs  being  optimized,  and  W is  the  weighting  factor  applied 


1 


1 
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to  each  output. 


It  should  be  noted  that  the  performance  measures  are  an 
approximation  and  when  the  weighting  factor  is  introduced 
the  performance  measure  diverges  significantly  from  the  defining 
integral  equation.  Comparison  of  performance  measures  for  a 
given  plant  with  various  types  of  compensation  should  be  done 
with  the  same  weighting  factors.  The  same  performance  measure 
could  be  achieved  for  a system  using  two  different  compensation 
schemes,  since  the  performance  measures  for  each  channel  may 
vary  significantly. 
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C.  A TWO  INPUT-T^TO  OUTPUT  SYSTEM  EXAMPLE  | 

1.  Transfer  Matrix  Solution 

In  order  to  demonstrate  the  program  it  was  decided  to  i 

analyze  a simple  system,  synthesize  the  compensation  required  | 

to  produce  a specified  output,  and  then  show  that  the  optimi-  - 

zation  program  will  arrive  at  the  same  results.  It  is  emphasized  ^ 

that  this  program  will  not  determine  the  type  of  compensation  j 

required  to  achieve  the  desired  response.  It  will  set  free 
parameters  in  the  compensators  (specified  by  the  user)  to  most 
nearly  produce  the  desired  response. 

i 

The  analysis  and  compensation  was  performed  using  the 
transfer  matrix  method  described  by  OGATA  (2) . The  uncompensated 
plant  and  the  compensated  plant  were  simulated  to  show  that  the 
compensation  achieved  the  design  rpecif ications . Certain  para- 
meters in  the  compensator  were  defined  as  free  parameters. 

These  free  parameters  were  offset  from  their  optimal  values 
and  constraints  were  introduced  which  restricted  the  range  of 
values  for  these  parameters  during  optimization . The  plant 
and  compensators  with  free  parameters  were  simulated  in  the 
optimization  program  and  optimization  was  accomplished. 

To  illustrate,  a simple  two-input  two-output  linear 
system  shown  in  Figure  3-1  was  selected  as  the  plant  to  be  , 

compensated.  The  design  specifications  chosen  were: 

1 . Channel  one  and  channel  two  are  to  be  decoupled  during 
the  transient  period  as  well  as  in  steady  state. 

2.  Channel  one  is  to  have  a second  order  response  to 

a step  input  with  a natural  frequency,  , of  10,  a 
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2: 


I 


damping  factor,  of  0.4,  and  no  steady  state  error. 


L. 


3.  Channel  two  is  to  have  a second  order  response  to  a 
step  input  with  a natural  frequency,  W^,  of  4,  a damping 
factor,  5,  of  0.6,  and  no  steady  state  error. 

To  achieve  total  decoupling  the  off-diagonal  elements  of 
the  transfer  matrix  of  the  compensated  system  must  be  zero. 
The  desired  compensated  system  transfer  matrix,  G(s),  is: 

I 100 0 

+ as  + 100 


G(s) 


0 


16 

+ 4 . as  + 16 


{52) 


The  input/output  equations  obtained  from  Figure  3-1  are: 


Y^(s) 


U^(s) 


(S+7) (S+12) 


^2(3) 


2 

(3+12) 


(53) 


Y2  (s) 


02(3) 


1 

(s+2) (s+4.3) 


+ Y^(s) 


4 

(3+2) 


(54) 


Substitution  of  equation  (53)  into  equation  (54),  and  equation 
(54)  into  equation  (53)  produces  equations  (55)  and  (56). 


(s) 


2U2(s) 


3Y^(s) 


^1^^^  “ (s+7)  (s+12)  ^ (s  + 2r(si4.3)  (s+12)"  ^ (s  + 2)  ('s+n) 


U2(s)  4U^(s)  8Y2(s) 

^2^®^  “ (s+2) (s+4.3)  ^ (s+12) (s+7) (s+2)  ^ (s+2) (s+12) 

Collecting  like  terms  in  equations  (55)  and  (56)  and  rearranging: 
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I 


Y^(s) 


s +14S+16 


U^{s) 


2U2(s) 


Ti+im 


^2(3) 


s +14S+16 


4U3_(s) 

(s+12)  (s+7r(s+2) 


U2(s) 


3+4.3) 


Multiplying  both  sides  of  equations  (57)  and  (58)  by 

(s+2) (s+12) 
s^+14s+16 

(s+2)U, (s)  2U, (s) 

Y (s)  = ^ + f 2 

(s+7) (s  +14S+16)  (s+4 .3) (s  +14S+16) 


Y2(s) 


4U^(s) 


(s+12)U2 (s) 


(s+7)  (s  +14S+16)  (s+4. 3)  (s  +14s+16: 


The  transfer  Matrix  for  the  uncompensated  plant  is; 


G^(s) 


(s+2) 2 ; 

(s+7) (s^+14s+16)  (s+4. 3) (s^+14s+16) 

4 (s+12) 

(s+7) (s^+14s+16)  (s+4. 3) (2^+14s+16) 


The  generalized  compfensation  to  be  used  is  shown  in  Figure 
3-2.  Equation  (41)  will  be  used  to  calculate  the  required 
compensation;  since  in  this  case  the  feedback  compensation 
transfer  matrix,  H(s),  is  the  identity  matrix. 


{ [9(s)"^  - i]Gp(s) 


Since  G^(s)  is  a 2x2  matrix  in  this  example  G^(s)  and  G(s)  will 
also  be  2x2.  By  inspection  (Equation  (52)),  G(s)  is  non-singular 


and  can  be  inverted. 
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k 

I 


I 


I 

1 


The  compensated  system  transfer  matrix  G(s)  is  given  in 
equation  (52).  The  inverse  of  G(s)  is  calculated  by  equation 
(62)  . 


G(s)‘^ 


COF 

G(s) 


(62) 


,T  , 


Where  COF  is  the  transpose  of  the  cofactor  matrix  of  G(s), 
and  |G(s)|  is  the  determinant  of  G(s). 

16  0 


COF 


s +4.8S+16 


100 


3 +8.-3+100 


(63) 


G ( s ) I = 


(16) (100) 

— « 

(s  +8S+100) (s  +4.8S+16) 


(64) 


G V 3 ) 


-1 


s i-8s  + 100 
100 


s''  + 4,8s+16 
16 


(65) 


^s2+8s 
100 


Cg(3)"^-i]  = 


s +4 . 8s 
16 


(66) 


The  plant  trcuisfer  matrix,  G^(s),  is  given  in  equation  (61)  as: 
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K 

1 


(s+2) 

T 


(s+7) (s  +14S+16)  (s+4.3)(s  +14s+16) 


(s+12) 


(s+7)  (s"‘+14s+16)  (s+4.3)  (s  +14s+16) 


(61) 


Multiplying  equation  (61)  by  equation  (66) ; 


[ G(s)'^-I]Gp(s) 


3(s-i-8)  (3H-2) 


2s  (3-*-8) 


100(3+7) (s^+14s+16)  100(3+4.3) (s^+14s+16) 


4 (s+4.8) 


s(s+4.8) (s+12) 


16 (s+7) (s^+14s+16)  16 (s+4.3) (s^+14s+16) 


(67) 


If  this  matrix  is  non-singular  it  can  be  inverted. 


_ 1 


[G(s)  *-I]G^(3) 


3 (3+8)  (s+2)  (s+4.8)  (5+12) 

(16) (100) (s+7) (s+4.3) (s^+14s+16)^ 


(68) 


3a  (s+4.8)  (s+3) 

(16)  (100)  (s+7)  (s+4.3)  (3^+14afl6)^ 


I [(3(3j"^  - I]Gp(s)  I = 


s^(s+4)  (s  + 8)  (s^+14s  + 24-8) 

(16) (100) (s+7) (s+4.3) (s^+14s+16)^ 

s^(s+4.3)  (s+8) 

(16) (100) (s+7) (s+4.3) (s^+14s+16) 


(69) 
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One  of  the  possible  realizations  of  this  compensation  is 
shown  in  Figure  3-3.  To  obtain  a confirmation  of  the  design, 
the  system  was  simulated  with  and  without  compensation  using 
International  Business  Machines'  Digital  Simulation  Language, 

DSL  (13) . Figures  3-4  through  3-9  show  the  results  of  that 
simulation,  which  indicate  that  the  compensators  used  do  in 
fact  achieve  the  stated  objectives  of  the  design.  In  each  plot, 
trace  number  1 is  the  desired  response,  Z,  and  trace  number  2 
the  system  response,  Y.  In  Figures  3-4B,  3-73,  3-8B,  and  3-9B 
the  two  traces  are  superimposed.  In  Figures  3-5A  and  3-6A 
the  desired  responses  are  zero.  The  reader  should  note  the 
scale  factors  present  on  some  of  these  plots. 


2.  Parameter  Optimization  Solution 

The  compensated  system  and  the  desired  response  curves 


were  then  simulated  and  solutions  were  obtained  using  the 


32 


Figure  3-6A 

Y1  and  Z1  vs  Time  with  Rl  = 0.0,  R2  = Unit  Step 

= 2.91  X 10"^) 
ss  ' 


- t 

'3.00  3.00  4.30  S.OO  1.00  10.00 

Figure  3-6B 

Y1  and  Z1  vs  Time  with  Rl  » 0.0,  R2  = Unit  Step 

For  Compensated  System 


Figure  3-8A 

Y1  and  Z1  vs  Time  with  Ri  = Unit  Step,  R2  = Unit  Step 


Figure  3-8B 

Y1  and  Zl  vs  Time  with  Rl  = Unit  Step,  R2  = Unit  Step 

For  Compensated  System 


Figure  3-9A 

Y2  and  Z2  vs  Time  with  R1  = Unit  Step,  R2  = Unit  Step 


Figure  3-9B 

Y2  and  Z2  vs  Time  with  Rl  = Unit  Step,  R2  = Unit  Step 
For  Compensated  System 


program  discussed  in  this  section.  The  time  response  results 
are  shown  in  Figure  3-10.  There  is  a slight  difference  between 
the  desired  response  curve  and  the  output  of  the  simulated 
system. 

Next,  two  parameter  optimization  was  performed.  The  pro- 
portional gain  constant  of  each  compensator  was  selected  as 
the  adjustable  parameters.  Each  of  the  four  performance  measures 
given  by  equations  (47)  through  (51)  were  used  in  successive 
optimizations.  All  other  aspects  of  the  program  were  unchanged 
during  this  optimization.  Equal  weighting  was  applied  to  each 
channel.  The  gains  100  and  16  were  chosen  because  their 
variation  should  have  similar  impact  on  the  two  respective  out- 
puts. The  upper  limits  were  105  and  21  and  the  lower  limits 
were  95  and  11  respectively.  The  integration  step  size  was 
.02  and  the  final  time  was  5 seconds.  Optimization  was  begun 
with  the  two  gains  set  at  105  and  11  respectively,  and  Rl  = R2  = 
Unit  Step  functions.  The  results  are  given  in  Table  3-1. 

The  system  was  simulated  using  these  values  for  the  two 
free  parameters.  The  simulation  results  are  shewn  in  Figures 
3-11  through  3-14.  This  simulation  indicates  the  time  multiplied 
integral  square-error  performance  measure  produced  the  best 
results . 
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Performance 

Measure 


Value  of 
Minimum 
Performance 
Measure 


^100 


G 


16 


Integral  , 

Square-  1.2026097  x lo"  104.9999  20.99998 

Error 


Time 

Multiplied 

Integral  1.641641  x 10~  102.2554  15.27091 

Square- 

Error 


Integral  , 

Absolute  1.353197  x 10*  104.9862  14.32275 

Error 


Time 

.Multiplied  . 

Integral  9.994698  x 10  103.2728  14.31198 

Absolute 

Error 


Table  3-1 


RESULTS  OF  OPTIMIZATION 


CC2  ooii  G06  008  OTO 

Figure  3-lOB 

Y2  and  Reference  2 vs  Tine  Simulated  with  Thesis  Program 
Parameters  Set  at  Analytical  Optimoim 
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Figure  3-llA 

Y1  and  Reference  1 vs  Time  After  Optimization  With 
Integral-Square  Error  Performance  Measure 
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Figure  3-113 

Y2  and  Reference  2 vs  Time  After  Optimization  With 
Integral  Square-Error  Performance  Measure 
(Curve  With  Greatest  Overshoot  is  ZZ) 
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Figure  3-12B 

Y2  and  Reference  2 vs  Time  After  Optimization  With 
Time  Multiplied  Integral  Square-Error  Performance  Measure 
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Figure  3-13A 

Yl  and  Reference  1 vs  Time  After  Optimization  With 
Integral  Absolute  Error  Performance  Measure 
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Figure  3-13B 

Y2  and  Reference  2 vs  Time  After  Optimization  With 
Integral  Absolute  Error  Performance  Measure 
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Figure  3-14A 

Yl  and  Reference  1 vs  Time  After  Optimization  With 
Time  Multiplied  Integral  Absolute  Error  Performance  Measure 
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CONCLUSIONS 


From  the  results  of  the  two-input  two-output  example 
system,  the  following  conclusions  can  be  stated: 

1.  The  technique  of  cost  function  minimization  in  the 
time  domain  can  be  used  effectively  for  compensator  parameter 
optimization  in  multiple  input  multiple  output  control  systems. 

2.  Cost  function  minimization  is  not  a substitute  for 
control  system  design.  The  user  must  be  able  to  select  the 
proper  type  of  compensator  which  is  capable  of  achieving  thej 
desired  results,  before  initiati.ng  the  optimization  of  the 
free  parameters . 

3 . The  four  performance  measures  evaluated  produced 
slightly  different  solutions  to  the  same  problem,  with  the 
ti.me  multiplied  integral  square-error  yielding  the  best 
results  in  the  example  given. 
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APPENDIX  A 


COMPENSATOR  □PTI'^IZATION  IN 
MULTIVARIABLE  CONTROL  SYSTEMS 
BY 


JOHN  T.  MOWREY 
MAJOR  USMC 
DECEMBER  1973 

THIS  PROGRAM  WAS  DEVELOPED  AS  A GENERAL  PURPOSE  DESIGN 
AID  TO  BE  USED  to  OPTIMIZE  FREE  PARAMETERS  IN  THE  COMP- 
ENSATCRS  OF  A MULTIVARIABLE  CONTROL  SYSTEM.  IT  CAN  ALSO 
BE  USED  FOR  THE  OPTIMIZATION  OF  FREE  PARAMETERS  IN  A 
SINGLE  INPUT  SINGLE  OUTPUT  SYSTEM'  OR  ^0  SIMULATE  EITHER 
TYPE  OF  SYSTEM.  OPTIMIZATION  IS  ACCOMPLISHED  BY  MINIMIZ- 
ING THE  ERROR  BETWEEN  A USER  SPECIFIED  REFERENCE  CURVE 
AND  THE  OUTPUT  OF  THE  SIMULATED  COMPENSATED  SYSTEM. 


DATA  CARDS 


CARD  1 NRUNSf  N2PT»NGRAPri,NREF,I  IN.NOLTtlPRINT  ,IFRE(3 
FORMATt  BIS) 

NRUNS-FLAG  for  number  OF  RUNS.  NRUNS  MUST  EQUAL 
1 WHEN  OPTIMIZING. 

NOPT-FLAG  FOR  OPT  IM  iZATION/SI  VULATIQN 

1 -OPT  I MI Z AT  I ON 
Z-SIMULATION  ONLY 

(IF  N0PT  = 1 THE  PLANT  WILL  BE  SIMULATED  AFTER 
OPTIMIZATION  AND  THE  REQUESTED  OUTPUT  PROCUCEC 

NGRAPH-FLAG  FOR  TY°E  OF  GRAPHICAL  OUTPUT  DESIRED 

1- PRODUCES  PRINTER  plQtS 

2- pRO  DUCES  VERSA  TEC  PLOTS 

3- NO  GRAPHS  PRQOUCED 

NPEF-FLAG  FOR  REFERENCE  CURVE 

1 - PRO  DUCES  CURVE 

2- NO  REFERENCE  CURVE  IS  PRODUCED 

IIN-thE  number  of  INPUTS  (UP  tq  25) 

NOUT-THE  NUMBER  OF  OUTPUTS  (UP  TO  3) 

jPRInjt-claG  for  tabulated  DATA  FOR  REFERENCE 
CURVE  AND  SYSTEM  pESPOnSE. 

1- PRODUCES  PRINTED  OUTPUT 

2- NO  OUTPUT  DA7A  PRINTED 

ISREQ-*^  REQUENCY  OF  PRINTED  OUTPUT.  IF  IFREQ=1D 
EVEPY  lOTH  DATA  POINT  ^ill  BE  PRINTED. 


CARD  2 T^7,DT,TF 

F0RMAT(3F10,5  > 
TA7-INITIAL  TIME 
DT-IN^EGR  ATION  STEP  SIZE 
TF-FINAL  time 
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CARO  3 NV  ,NAV  ,r>rrA,NPR,  I ? , I PM 
F-J^mak  sis  ) 

this  card  must  3E  CMITTED  IF  KaP'^=2 

NV-THE  NUMBER  OF  FREE  VARIABLES  '^3  BE  SET  BY 
HPTIM  IZATIDN. 

NAV-THE  NUMBER  OF  AUXILLIARY  VARIABLES 

NTA-THE  NUMBER  OF  TRIALS  ALLOWED  BOXPLX 

NPR-FREQUENCY  OF  OUTPUT  FROM  BOXPLX  FOR 
DIAGNCSTIC  PURPOSES 

IP-FLAG  FOR  I NTEGER /REALMS  OPTIMISATION. 

* « 

* NOTE  * 

* * 

SEE  BOXPLX  FOR  ADDITIONAL  INFORMATION  ON  THE  5 VARIABLES 
LISTED  ABOVE.  BOXPLX  WAS  CONVERTED  TO  DOUBLE  PRECISION 
= OR  USE  IN  THIS  PROGRAM;  HOWEVER  the  COCU  me  NT  AT  I ON  WAS 
NOT  changed 

IPM-clag  for  PERFORMANCE  MEASURE. 

the  performance  McASiJPES  US  50  ARE  PROPORTIONAL 
TO  THOSE  INDICATED  BELOW. 

1- IMTEGRAL  square-error  (WEIGHTED) 

2- TIME  MULTIPLIED  INTEGRAL  SQUARE-ERROR 
(WE  IGHT  ED  ) 

3- INT=gral  absolute  ERROR  (WEIGHTED) 

A.-TIME  multiplied  INTEGRAL  ABSOLUTE  ERROR 

(WEIGHTED  ) 


CAR  DS  A— T XS  ( I ) 

FORM-A'iaFlO.S) 

ALL  CF  these  CARDS  ^ijST  BE  OMITTED  1=  NOPT=2 
XSd)  ARE  the  starting  VALUES  OF  THE  FREE  =APA- 
MFTP9S.  the  number  of  CAROS  REQUIRED  IS  DEPEN- 
DENT ON  NV.  ONE  VALLE  MUST  BE  READ  IN  FOR  EACH 
FREE  PARAMETER.  IE  IF  NV= B ONLY  CARO  A IS  RE- 
QUIRED AND  3 PARAMETER  VALUES  SHOULD  BE  ON  IT 
I N 8 = 10.  5 F OR  MAT. 


CAROS  8-11  XU(I) 

FORMAT  (8F10.5) 

these  cards  are  THE  SAME  AS  CARDS  A-7  EXCEPT 
XUd)  ARE  THE  UPPER  LIMITS  ON  THE  FREE  PARA- 
METERS. OMI T IF  N0PT=2. 


: ARC  S 12-15  XL(  I ) 

F0RMAT(3F10.5  ) 

these  cards  are  the  same  as  cards  4.-7  EXCEP’ 
XL(I)  ARE  THE  LOWER  LIMITS  ON  THE  FREE  PARA- 
METERS. OMIT  IF  NCPT  *2. 
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CARDS  16»18t  AND  20  IDUT, I WT, I ta3 

= DRMAT(  3F10.5) 

IQU'^-THS  NUNBER  OF  THE  BLOCK  FROM  WHICH  Thr  OUT- 
PUT 13  TAKEN,  I F 2 OR  MORE  BLOCKS  FEED  AN  OUT- 
PUT NOOEf A TYPE  1 BLOCK  WITH  G=1  MUST  BE 
PLACED  BETWEEN  THAI  NODE  AND  THE  OUTPUT. 

FOR  A 3 OUTPUT  SYSTEM  ALL  3 CARDS  ARE  REDLIRED. 

FOR  A 2 OUTPUT  SYSTEM  15  AND  18  ARE  REQUIRED. 

IWT-AN  INTEGER  WEIGHTING  FUNCTION  WHICH  WILL 
BE  CONVERTED  TO  R6AL#8.AMD  WHICH  ALLOWS  THE 
USER  TO  PENALIZE  THE  ERROR  BETWEEN  THE  SYSTEM 
OUTPUT  A.N0  THE  REFERENCE  RESPONSE  MORE  HEAVILY 
AT  ONE  OUTPUT  THAN  ANOTHER.  INTEGERS  BETWEEN  I 
AND  10  ARE  RECOMMENDED.  WEIGHTING  SHOULD  BE 
REDUCED  IF  OVERFLOW  IS  ENCOUNTERED. 

ItaB-A  FLAG  FOR  TABLE  LOOK  UO  OF  the  REFERENCE 
C URVFS. 

1- OROGRAM  WILL  READ  TABULATED  DATA  FOR  -he 
REFERENCE  CURVES. 

2- PROGRAM  WILL  CALCULA'^E  A SFCCNC  ORDER  RES- 
PONSE CURVE. 


CAODS  17,19»  AND  21-B  ET  A » DELTA  , WN , AMO  , DEL  AY  , I NP  LT 
cqRMATI  5P10.5  ,aA.) 

'■HESS  CARDS  INPUT  THE  DATA  RECUIREQ  TO  COMol''^E  A 
SECOND  ORDER  RESPONSE  CURVE.  1=  THE  USER  DESIRES 
THE  TABLE  LOOK  JP  FEATURE  THESE  CAROS  ARE  RE- 
PLACED BY  TABULATED  DATA.  IF  A TABULATED  REF- 
ERENCE CURVE  IS  DESIRED  FOR  OL'^PUT  NC.  1,  CARO 
17  IS  REOLACEO  BY  TABULATED  DATA  IN  SFIO.  5 FORMA 
•^HE  number  of  data  POINTS  REQUIRED  IS  T F/ Dj  (TF 
AND  OT  FROM  CARO  2)  . 

Cd)  BETA!  I) 

Rd)  ~ ( S**2>2*0EL“A(  I ) *WN(  I)»Si-wN  ( I )#’*2) 

C ri  ) IS  THE  DESIRED  REFERENCE  RESPONSE  ASSOC- 
IATED WITH  OUTPUT  lOUTd). 

Rd)  IS  THE  FORCING  FUNCTION  SPECIFIED  BY: 
AMO(I)-THE  AMPLITUDE  OF  THE  FORCING  FUNCTION 
DEL  AY  ( I )-OEL  AY  A ER  T 47  BEFORE  THE  FORCING 
= UNCTION  IS  AFPLIED(TIM£  DOMAIN). 
riPUTd)-THE  TYPE  OF  FORCING  F UNCTION  : S T£0  , 
RAMP,  DR  RARA. 

3ETA(I)-THE  GAIN  OF  THE  tpaNSFFR  FUNCiTON. 

delta(I)-the  damping  factor. 

WN(I)-TH£  NATURAL  FREQUENCY. 


CARO  22  N 

= ORMAT(  12) 

THE  NUMBER  OF  BLOCKS  IN  THE  SYSTEM  (UP  TO  25). 
CAROS  23-47  IC(  I ) , I0(  I ),  I£(  I ),  IVd  ) ,G  (I  ) fR  ( I ) tZ(  I ) 

FORmaT(  4I2f  2X,2F1C.5) 

the  number  DF  CAROS  REQUIRED  IS  DETERMINED  BY  N. 
IC(I)-th£  NUMBER  OF  THE  3L0CK(l-25). 
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ID(I)-THE  TYOE  3F  3L3CK.  SIX  TYPES  ARE  AVAILfiBLE 

IE(  I ) ivu ) 

iD(n=i  1 3 I 

^ * TII^  I IV  ( I ) 

1 S+P  I 

0(  ) - j S+Z  j IVU) 

I S+P  r ~ 

IE( I ) I 3 I I V( I ) 

I D ( I ) =4 

S**2+2*DELTA*V.M=>‘S  + WN**2 


ID( I ) =5 


I V(  1) 


I7(  I)  =3*  IE(  I ) FOR  O^lEd)  LE  O. 

AMO  0*IE{  I ) GE  LL 
17(1)  =UL  FOR  G^IEd  ) GT  LL 

I7(  I ) = LL  =0R  3*IE(  I ) LT  UL 


I0.(  n =6 


5 / 


'I 


17(1) 


I£(  I ) 


/ -3 
/ I 3 
I/-I 


I 7(  I ) =G«IE  ( I ) =GR  G*I£(I)  GE  6 
OR  LE  -3 

17(11=0  FOR  3*IE(n  L"^  3 AMO  GT  -3 
IE(I)-THE  IMPUT  mode  MUM9E^. 

I7(I)-THE  OJTP'JT  MODE  MU'^BER. 

G(n-THE  GAIM  FIELD.  3(1)  REPRESEMTS  THE  GAIN 
OF  THE  BLOCK. 


P(  I )-thE  POL  E 
FOR  I0( I ) =2 
F^R  I 0( I ) =4 
FOR  ID(I)=5 
= 0R  I0(  1 1 =6 


FIELD  . 

OP  3 P(I)  = R0LF 

p( I )=dampimg  factor 

P(  I ) = UP'RER  L IM  IT,UL. 
P( I ) ^REFERENCE, 8. 


Z( I )-th£  zero 
FOR  10(1 ) =3 
FOR  I0(I)=4 
FOR  I 0 (I)  «5 


FIELD 
Z( I ) *ZER0 

Z(  I l = MAT'JP  AL  FR  EQIJEMC  Y ,WM. 
Z(  I ) SLOWER  LI^^I  T,  LL. 
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CA50S  4a-r72  AMP  (I  ) , 0ELAY(  I ) , INPUT(  I ) , I5E(  I I 
PQRMAT (2F10.5 ,AA»  12 ) 

ONE  CARO  IS  NEEDED  FOR  EACH  FORCING  FLMCTION  ON 
THE  SIMULATED  COMPENSATED  SYSTEM. 

AMP(  I),  DELAYm»  ANDINPUT(  I ) APE  AS  DESCRIBED  FOR 
CARDS  17, 19, AND  21, EXCEPT  THEY  APOLY  TC  THE  COM- 
PENSATED SYSTEM  IN  THIS  CASE. 

IEE(I)-TH5  NODE  TC  XHICH  ^hE  FORCING  FUNCTION  IS 
APPLIED. 

CARDS  73-78  TITLE  CARDS  FOR  VERSATEC  PLOTS.  IF  NGRAPH=2, 
NOUT+2  TITLE  CARDS  ARE  REQUIRED(2  PER  OUTPUT) 
TITLES  GO  IN  COLUMNS  1-48  ONLY. 

« * 

♦ NOTE  * 

■m  ♦ 

WHEN  CPTI  MI  ZING  ,ONc  FORTRAN  CARD  OER  FREE  VARIABLE  MUST 
3E  INSERTED  IN  SUBROUTINE  PLANT.  THE  APOROPRIATE  POSITION 
IS  IDEN-^IFIED  BY  THE  COMMENT  CARD  » ENTER  G(I)=C(I)  VAL- 
UES AT  TH I S POI  NT  ' . 


The  follow  ING  JCL  cards  should  follow  JOB  CARO,  EACH 
BEGINNI  N3  I N COLUMN  1, 

//  EXEC  FCR'^CLGV,REGI0N.Ga  = 350K 
//FOP  T.  SYSPR  I NT  DO  DUMMY 
//SYSIN  DO  * 

THE  SECOND  OF  THESE  SHOULD  EE  REMOVED  IF  A PROGRAM 
LISTING  IS  DESIRED  AS  WcLL  AS  OUTPU*^  DATA. 


I'^PLICI'^  3EAL*8  (A-H,0-Z) 

REAL‘3  XS(25),XL{25  ),XU(25),X(  2),XC0T(  2),DELTA(  3), 
*W(3),«N(3)  , beta (3) ,THA CUT (3001 ,3)  ,XOATA ( 300  I, 3 ) , 

1A'^0(3  ),0£LAY(3) 

DI MENS  ION  IDUT(3 ) ,IWT (3  ), ITAB (3) , INOUT( 3 > 

INTEGER  S'^E3,RAMP  ,?aRA 

DA'^A  ST  cO,R  AMP  , PARA /‘♦HSTEO,  A-HP  AMP  ,4HPARA/ 

COMMON  THAa'JT,XOATA.T,OT,TF, 

*NOUT,  I IN  ,NV  ,M3,IC0NT,  NEC,I  SKI  P , I TF,iauT  ,NCPT 
COMMON  /REGl/  W,  IPM 

INPUT  CONTROL  FLAGS 

RE  AO  ( 5, 5 0 NRUNS, N OPT,  N GRAPH, NR  SF,  1 1 N , NQUT , I PR  INT  , 1 FR  EO 
WRITE!  6,51  )NRUNS,  II  N,NOUT 
IF(NGPT.E0.1)WR  IT  £(  6,52  ) 

I =(  NOPT.eq.2)  wRI  TE(  6 ,53  ) 

I F(NGRAPH.5Q  .1  )WR  ITE(  6, 54) 

IF(NGRAPH.E3.2)  WRITE(6,55) 

IF(N3RAPH.Eq.3) WR ITE( 6,56) 

IF(NREF.EQ.l  ) ^R  IT  £(  6,  37  ) 

IF{NREF.Eq.5)wRITE(6,5a) 

IFdPRINT.EO.l  )WRITE(6,60)  IFRE3 
IF(  IPRINT.EC.2  )WRITE(6,59) 

''O  30  IP  UN=1  ,NRUNS 
’EAO(  5,61  )T47,DT,  TF 
itf»tf/dt 

IF(  ITF.GT.3001)  I TF  = 3C01 

IDP  s ITF 

NPPLT*I0P/5 
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33 

32 


35 
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WRITS(6,62)TA.7,0T,TF 

INPUT  QPTINIZATICN  DATA 

GD  TD  (1,2),N0PT 
READ  (5 *50  »NV,NAV,NTA,NPR, IP,  I PM 
REAOC  5,61  )(  XS(  I I ,1=1  ,NV) 

REA0(5,61  )(  XU(  I >,I=1,NV) 

REA0(5,61  XXL  (I  ),I=!,NV) 

WRI TE( 6,6  3) NV,NA V ,NTA ,NPR,I  P ,I PN 
DO  31  1=1, NV 

WRIT5(S  ,64) I,XS(I ) ,I ,XU (I ) , I ,XL ( I) 

CONTINUE 

GO  TO  (3,4),NREF 
00  32  I=1,N0UT 

input/calculate  REFERENCE  CURVES 

READ (5, 50>I0UT( I) ,IWT( I) ,ITAB<I ) 

W(  I ) = IWT ( I ) /I  00000. CO 
WRI TE(6 ,65)  I ,iaUT(  I ) ,I , W ( I ) 

IF(  I-^A8(  I)  .EQ  .1  IWRITE(  6,76)  I 
IF(ITAe<n.E0.2)WRIT£{6,77  )I 
ITA=ITA8<I ) 

GO  TO  (5,6),I'^A 

REAC(5 ,61  ) (XDATA(  IDATA,  I)  , IDATA=1,  ITF) 

GO  TO  32 

0 EA0(5,?5  )9ETA(  I ) ,0ELTA(  I)  ,WN(  I ),AMP  ( n ,0ELAY(  I ) , 
INPUT! I ) 

WR  ITE  ( 6,66)  I ,3ETA  ( I ) , I,DE  L’A  ( I ) ,I  , <.N  (I  ) ,I  ,AMP(  I ) , 
1,  INPtJT  ( I ),  I,  DELAY  ( I I 
T = T47 
A17=0.D0 

1 DAT  A=0 
X(1  ) = 0.00 
X ( 2 )=0.D0 
NT»0 

03  33  J *1  I TF 

IF(".Ge!d£LAY(  I ).  AND.  INPUT!  I)  . EQ  . STEP  ) A1  TsAMP  ( I) 
I !=(T.GE.OELAY(I  ) .AND.  IN  OUT  ( I ) . EQ.RAMP  )A17=AMO(  i ) 
’-DELAY!  I ) ) 

IF!T  .GE.DELAY!  I ) .AND.  IN  PUT!  I ) .EQ.PARA  )A17  = AMO(  I) 
» ( (T-DELAY!  I)  )**2 ) 

XOOTI 1)=X!2 ) 

XD.0'^!2)=-2.C0»C£LTA!  I)*WN!  I )*X!2)-!WNII  )*=2)* 

X!1  X-3ETA!  I )=A17 
S = RKL0EQ!  2;  X,  XDOT,  T,DT,NT) 

IF!S-1. 00)10,11, 16 
WR  ITE!  6,6  7)  I 
STOP 

XOATA! J,I ) =X!1 ) 

continue 

CONTI  NUE 
GO  TO  12 
DO  34  1=1, ITF 
DO  35  J*1,N0UT 
XDATA!  I , J)  = O.DO 
CONTINUE 
C ONT  I NU  E 
T = T4T 
ISKIO  = 0 

C 

l<=  SIMULATION  ONLY  CALL  PLANT  C 

C 

IF!  NOPT.EC.2  ) CALL  PLANT!C) 

IF!NaPT.EQ.2)  GO  TO  13 

C 

optimization  C 

c 

o«l  .DO/3  .00 
WRITE  <6,73  ) 

CALL  BOX?  LX!  NV,  NAV,  NPR,NTA,  R,  XS  ,I?,XU  ,XL,YMN,  I£,R  ) 
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WRITS (6t73) 

WRITE(6f63) 

30  36  I»1»NV 

WRITE! 6f  69)1,  XS( I ) 

36  CON-^INUE 
WRITE(6,70)  YMNflER 

PRINTEC  OUTPUT 

3 GO  70  (14,15),IPRINT 

4 WRITE(6,'^iT 
00  37  Ial,ITF ,IFREQ 

T P*0T*I 

WRI TE (6 ,72)  TP, (XOATAI I, J) , Th ACUT ( I , J ) , , NOUT ) 

37  CONTINUE 

PLOTTED  OUTPUT 

15  GO  TO  (7,9,9)  ,NGRAPH 

7 00  38  II=»l,NOUT 

WP  ITE(6,73) 

CALL  POLT  (NPPLT  , lOP,  II  ) 

33  CONTINUE 
GO  -0  9 

a 00  3 9 II  =1  ,NOUT 
WR  I TE(  6,73) 

CALL  PIC(NPPLT, lOP,  I I) 

39  CONTINUE 
9 HRITE(6,74)  IRUN 

30  CONTINUE 
STOP 

50  FORMAT! 815) 

51  FCRMATCOS  I2,1X,'R'JN(5  ) WILL  RE  ^A  OE  ' , 4X  , I 2,  IX , 

* ' INPUTS’ ,4X,  I 2, IX,'  OUTPUTS'  ) 

52  forma-! '0*  , 'OPriM  IZAT  ION  CALLED  FOR') 

53  FCRMAT!<0’  , 'SIMULATION  CNLY  CALLSC  FOR') 

54  FORMAT!  'O' , 'PRINTER  PLGTIS)  CALLEC  FOR') 

55  FCRMA-^  ! *0*  , 'VE.RSAT  EC  PLOTS  CALLED  FOR') 

56  =0RMAT(  • O'  , ' NO  GRAPHICAL  OUTpU”  C^LLEO  FOR') 

57  FOBM  at  ( *0',  'REFER  EMC  E CURVE  CALLED  FOR') 

53  FORMAT  <' O' ,' NO  REFERENCE  CJO'/E  CALLED  FOR') 

59  poRMAT!  'O'  , 'NO  PRINTED  SIMULATION  DATA  CALLED  FOR') 

60  FORMAT!  'O' , 'SIMULATION  DATA  PRINT  FREQUENC Y='  , 1 X , I 2 ) 

61  FORMAT!  8F10. 5 ) 
bt  FORMAT!  'O' ,' INI  TI  AL  T I M E*'  , FI  0 . 5 , 3 X , ' STE®  S I Z E=» ' ,F  1 3 . 5 

«,3X,'FIN-AL  T IME=' , F10.5  ) 

63  = CRMA''’ ! 'O'  , 'NV^'  , 1 2 , 3X  , ' NA  V =’  , 12  ,3  X , ' NT  A=  ' , 15 , 3X  , ' NPR  a 
*',I5,3X,  'IPa',Il,  3X,  ' IPM  = ' ,11) 

64  format  ! 'O'  , 'XS(  ',  12, ' ) =*  ,F1  0 .5,  5 X,  ' XU!  • , I 2,  ' ) =»  ' , FIO.  5 , 
»5X,'XL(  ’ ,12,'  )a'  , FI 0.5) 

65  FORMA-!  'O’  , ' lO'JT!  • ,11,  • )a'  , 12, 5X,  'W!  ',!  1,  ' ) = ' ,F10.  5) 

66  FORMAT! 'O'  , ' BET  A!  ’ , II , ' ) »' , =10 . 5 , 3X , ’ DEL T A ( • , II,  ' ) = ' , 
♦F 10. 5,3X, ' WN!  ’ ,11  ,'  ) a'  ,F10. 5,3  X, ’ AMP! • , 11  ,'  )»•  ,Fia  .5, 
*3X,'  INPU'''!  ',  II,  • )a',  AA  , 3X,  'DELAY!  ' , I 1,'  )a'  ,F1C.5) 

67  fqrma-! 'O'  ,' INTEGRATION  T ROJ  3L  S-R  EFSR  EN  C = CURVE') 

68  FORMAT!  '0 ', 'THE  FREE  PARAMETERS  ASSOCIATED  WITH  THE 

• minimum  PERFORMANCE  MEASURE  ARE:') 

69  F CRmaT(  'O'  , 'XS!  ' ,12  ,'  ) «'  ,F20.1  0) 

70  F^'SMiT  ( 'O' , 'THE  MINIMUM  PERFORMANCE  MEASURE  IS:',2X, 
•F20,10, 5X, ' I EPa'  ,2x,  II  ) 

71  FORMAT!  *1'  ,3X,’  TI  ME'  , 9 X , ’ XD  A- A ! 1 ) ' , 7X  , ' TH  AQU- ( i ) ' , 

* 6X,  ' XOATA!  2) ' , 7X,  ' THACU-!  2)  ' , 6 X , ' XDA  ta  ! 3)  ' , 7X,  ' 
*ThaOUT(3  )'  ) 

72  FORMAT!  'O'  ,6!  FI  0.5,  3X), FI  0.5) 

73  forma-!'!') 

7-^  =CRMAT(  ’ O'  ,' RUN  NUMBER'  , lx  , 11  ) 

75  FORM  AT  ( 5F10  .5,A4) 

76  format  ('O' ,' REF. CURVE!' , II,  •)*»,  '■'■ARULATED  DATA') 

77  FORMAT!  'O'  , 'R  £F  . C JR  YE  ! ' , 1 1 , ' ) * S EC  CNO  ORDER  RESPONSE') 
END 
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SUBROUTINE  PLANT (C) 

SUBROUTINE  PLANT (C)  SIMULATES  THE  SYSTEM 
implicit  REALT3  {A-H,a-Z) 

REALMS  G(25)  ,P(  25  » ,Z  ( 25),0RI  VE  ( 25),  THAI  25)  ,OMG(  25)  , 
ITHA0QT(23 ) , CMGOQT (25 ) , CRV I N ( 25 ) , = LAG! 25 , 25 ) , 

2 X2(25),  X2DOT(25),  THACUTI 3 001  ,3  ) , 

3XOATA(300l, 3 ), C(25)  , AMO( 25 ) , DE L A Y ( 2 5 ) , TAG ( 25,  25) 
DIMENSION  IC(25),ID(25),I5(25)  ,MF(25)  ,IR(25),  IV(25), 
IIEE(25),INPUT(25)  ,I0UT(3) 

INTEGER  STEP,  RAMP,  PARA 

DATA  STEP,RAMP,PARA/^hSTEP,AHRAMP,AHPARA/ 
common  THA0UT,X0ATA, T,DT,TF, 

CNCUT,  II  N,NV,M3,  ICONT  ,NEQ,  ISKI®,  ITF,  iaUT,NOPT 
IF  ( I SKI  P-1  ) 1,5,5 
I I SKIP  = 2 
REA0(5,2^)  N 
HI  » OT 
h2  =»  0.500*H1 
Nil  = 0 
N55  = 0 
M66  = 0 
ICK  = 2*N 
WRITE  (6,310) 

00  3 1 = 1, N 

INPUT  BLOCK  DATA 

READ ( 5,25  ) IC  ( I ) , 1 0 ( I ) , IE ( I ) , IV ( I ),G( I ),  P( I ) . Z(  I) 
IF(  I0(  I).  EO.l)  NilsiNll^-l 
IF(  I0(  I ).E0.5  ) N55=N55^1 
IF( ID( I) . EQ.6)  N66=N66>1 

WRITE ( 6, 26)  IC( I ) , !□( I ) , IE( I ) , IV( I ) ,G( I ) ,P( I ) ,Z( I ) 
3 : CNT  I MJ  £ 

WRI  TE  (6,312) 

00  316  1 = 1, IIN 

INPUT  FORCING  ^UMCTICNS 

READ!  5,311)A.M3(  I ) , 0£LAY(  I ) , INP'J*^(  I ) , lEE  ( I ) 

■=RITE(  6,317)  IE£(  I)  ,AMO(i)  ,implT(I)  , DELAY  ( I ) 

6 C'^NTINUE 

NJll  3 4«Nll 
n;5  = 4=n55 
N66  = A.MN66 
NEQ  = N-1 

SET  POINTERS 

DC  A J*1,N 
00  ^ K«l, N 

= LAG  (K,  J)  *0.00 

IF(  lV(K).eQ.IE(  J)  )FLAG(K,J)=1.00 

A continue 

INIT  lALIZATI  ON 

5 00  6 ICLR=1  ,M 

THAI  ICLP  ) = 0.00 
THA00T(ICLP)=0.00 
OMG( ICLR)=0. DO 
OMGOOT  ( ICLR  )*0.00 
X2  ( I C L R ) = 0.  00 
X200T(  ICLR)  =0.0  0 

6 continue 

T*O.CO 
NR2  * 1 
NR3  * I 
NRA  a 1 
M 2 » 0 

ICONT  3 0 
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IWAIT  = 0 

I T » 0 
ILAST  = 0 
I5LAST  = 0 
IfcLiST  = 0 

IFCNOPT  .EQ  .2)  50  TO  7 

ENTER  C(I)=G(I)  VALUES  AT  THIS  PCIN-^ 

7 DC  9 MDRV=1,N 

DRIVE(MDRV) =0.00 
ORVIN  (NDRV)=0.00 
0 0 8 M=1  , N 

IF( IV(M).NE.IE(MORV) ) GO  TO  SJS 

DRIVE  (MDRV»  = 0R1V  E(MCRV  )+THA(m1*FLAG(  N,MOR  V» 

315  I^dSEINKNE.  lE(MOPV))  GO  TO  P 

IF {INPUT(M) .EQ.STEP)  go  TO  32 
IF( INPUT (MJ .£Q. RAMP)  GOTO  33 
I=(INPUT(MI  .EQ.PARA)  GO  TO  3A- 

32  IF( T .GE.OELA Y(M I )DRVIN{M0RV ) aAMP( M) 

GO  TO  8 

33  IF(T.GE.OELAY(M» ) ORVIN (MQRV I aAMP( M»* ( T-OELAY (M ) ) 
GO  TP  3 

34  I F (T.GE.OELAY  ( *<)  ) ORV  in  (MORV  )=  AMP(m  ) *(  (T-0ELAY(  m ) ) 

♦ **2) 

3 CONTIM'JE 

0RIVE(»<0RV)  *0RIVE(  '*ORV  ) +ORV  I M (m  CRV ) 

9 CONTINUE 

10  M3  * M3fl 

SLOCK  SIMULATION 

I = (IWAIT.  EQ.O)  T = T+H2 
IF  (I  0(M3  ) .£0 .1 ) GO  TO  U 

IF  ( I0(  M3  ) ,£0^  I GO  TO  12 

IF  ( I0{ M3I.Ea.3)  GO  TG  13 
IF  ( ID( M3  ) .EQ  .4)  GO  TO  14 

IF  (I0(  M3»  .EQ.5  ) GO  TO  15 

IF  ( I0( M3) .£Q.6>  GO  TO  U 
WRI-E  (6t23) 

STOP 

11  THA(M3)  = G(  M3)*0RIV£(  M3) 

IWAIT  a IWAIT  1-1 

ILAST  = ILAST>1 

IF  ((  ILAST. £Q. Nil)  .AND.  ( ICONT. EC. NEw)  ) GO  TO  21 
GC  ^0  13 

12  THA00T(M3)  a -P  (M3)*ThA  (M3  ) •►G(  M3  ) *C?IVE  ( ><3  ) 

3 = RKLD£2(  THA,  THADGT.NRZ) 

IWAIT  a iWAITfi 

I F(  S-1. 00)1  7,13, 23, 

13  0MG00T(M3)  = -P  (M^)=0MG(  M3  ) -t-GI  M3)  *0RI  VE  (M3) 

S = RKL0E2  ( CMG,  OMGOOT  ,NR3) 

THA(m3)  3 ( Z(,M3)-P(  M3)  )wGMG(M3  )>G(M3)*0RIV£(M3  ) 

IWAIT  a rwAIT+i 
IFIS-l. 00)17,13, 21 

14  V a ::plX(  P,Z,G, DRIVE, X2,0MG,NR4) 

ThA( M3  ) s 0MG(M3 ) 

I WAIT  a IWAl'^4-1 
IF(V-1,0C)  17,  13,21 

15  THA(m3)  a G(M3  )MORIV  E(  M3  ) 

1=  (THAI  M3  ) .GT.  P(  M3  ) ) th  A ( M3  ) aP  ( m3  ) 

IF  (’HA  (M3  ) .LT  .Z(  M3  ) ) THA  ( M3  ) = Z ( M 3 ) 

IWAIT  a IWAITH 

I 3LAST  a I 5LAST+1 

IF  ( ( I5LAS’.£g.N33)  .AND.  (ICONT.EQ  .NEO)  ) GO  t^  21 
GO  TO  18 

16  THA(m3T  a G(M3)*0RI V£( M3) 

IF((0A3S(THA(M3))  ).LT.P(M3)  ) THA  ( M 3 ) aQ  .0  C 
IWAIT  a IWAIT<-1 

IoLAST  a I6LAST4-1 

IF  ((  I6LAST.EQ.N66).ANC.(  ICONT.EQ.NEQ))  GO  TO  21 
GO  TO  13 
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17  WRI  T=  (6,29) 

STOP 

13  ICuNT  a ICGNT>1 

1=  (ICONT-N)  10,19,20 

19  •')R2  * NR2+1 
NR3  * NR3+1 

3 NRA^I 
M 3 a 0 
I CONT  a 0 

19  ( IWAIT.EO.ICiO  I vi)AIT«0 
GO  TO  7 

20  v^RITE  ( 6,30  ) 

STOP 

STORE  OUTPUT  DATA 


21 

IT  a IT+1 

DC  33  lal 

TH  A OUTT 

38 

CUNT INUE 

IF  (T-TF) 

22 

NR2  a 1 

T,n  aTHAdOUTCI)  ) 

22,22,23 

M93  = 1 

NR4  a 1 
*<3  a 0 

ICON-'  a 0 

I WAIT  3 0 
ILAST  a 0 
15LAS-'  a 0 
I 6LAST  a 0 
GO  TO  7 

23  RETURN 

24  FORXATt  12) 

25  forma-' ( 412,  2X,3F1  0.5  ) 

26  =CRMAT  (/,5H  5Li^  ,i2,Il,2H 

27  format  (//,  2X,  • theta  OUT  IS  tha(',I2,') 

23  FC9MA-'  (40H  EON  SWITCH  CONTROL  DID 

29  FORMAT  (20H  INTEGRATION  TROUBLE) 

30  FOR'^ATt  sax, ’ERROR  IN  INTERGA’'!  ON.  • , 2X, 
I'ATTEMPTED  TO  INTSG.  '<CRE  THAN  N-EC’N'  ) 

3 10  forma* ( 'O' , 40X,  ■ ' 

311  format  (2F10.5,  A4,  12  ) 

312  FORMAT! ' C ,40X, ' THE  SYSTEM 
317  forma*  * " 

E NO 


,212,6X,3E20.7) 

~ •.rj.iji) 


NOT  WORK  **a) 


) 


FORCING  FUNCTIONS  ARE:') 
( ’O’,  'CR7IN(  ',  12,  • )a’,F10.5,A4,F10.5) 


FUNCTI ON  RKLDEQIN ,X,XDOT,T, OT,NT  ) 

CALCULATES  RESPGNSc  OF  SECOND 
ORDER  REFERENCE  CURVE 

ImoliOIT  REAL«8  (A-H,0-2) 

REAL*0  X(2),XXT(2),3(25) 

NT  a NT  ♦ 1 

GO  TO  ( 1,2,  3,4)  ,NT 

1 Hi  a DT 

H2  a Hi*C.500 
h3  a Hi  *2.300 
H6  a H1/6.0DO 
DC  11  J a i,N 
11  0(J)  a o.DO 
A a 0,5  00 
T a r ♦ H2 
GO  TO  5 

2 A a 0.  P923932139134525 
GO  TO  ? 

3 A a 1 .7071067311365475 
T a T * H2 

GO  TO  5 

4 DC  41  I a 1,N 

41  X(I)  a X(I)  ♦ H6-»X00T{I)  - Q(I)/3.00 
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NT  a 0 
RKLDEO  = 2. 

GC  TO  6 

5 DC  51  L = ! »N 

X(L)  = X(L)  ♦ 0T*XD0T(L)-0(L)  ) 

51  Q(L)  = H3*A*XOOT(L)  ^ (l.CO  - 3.00*A)aQ(L) 
RKL06Q  a 1. 

6 RETURN 
END 


FUNCTION  RKLDE2  (X,XD0T,NR2) 
implicit  REAL*3  (A-H,a-Z> 

REALMS  THAQUT(3001,3)  ,XDATA(  3001,2) 

REAL*?  X(<V),  XDOTI^I,  C(25) 

0  IMENSION  IOUT(  3) 

COMMON  THAOUT,XOATA,T,DT,  TF, 

CNCUT,II  N,NV  ,M3  tlCQNT  , NEC,  ISKI3,  IT  F,  lOUT 
GO  TO  ( 1,2, 3,4) , NR2 

1 HI  a 0^ 

H2  a H1*0.5D0 
H3  a h1*2.000 
H6  a H1/&.000 
Q(M2)  a 0.00 
A a 0 . 50  0 
GG  rn  5 

2 A a 0.2^28932138134525 
GO  TO  5 

3 A a 1 .7071067811865475 
GO  TO  5 

4 X(M3)  s X(  M3)-t-H6»X00T(M2)-Q(M3) /3.0c 
RKL0E2  a 1. 

IF  ( ICONT.EG.NEQ)  Ri<L0E2a2. 

GO  TO  6 

5 X(M3)  * X(M3)  ^■Aa(  CT^XDCT  (M3)-Q{ '^3  ) I 
Q(M3)  a H3*'A*X00T(  M3)+(1.00-3.00*A)aO<  M3) 
PKL0E2  a 1, 

6 R:TURN 
END 


FUNCTION  R:<L0E3  (X,X0GT,NR3) 

IMPLIC  IT  0EAL»3  (A-H.O-Z) 

REAL*  8 THAQUT  (3301  ,3  ),XCATA(3001,2) 

REAL*9  X(4)  , XOOT( 4)  , Q(  25) 

OI^ENS  ICN  ICUT(3I 

CCMMCN  THAOUT,XOATA,T  ,CT  ,^=  , 

CNOUT,  I iN,.NV  ,M3,  IC  ON  7 , NEQ  , I SX I o , ITF,  ICUT 
GC  tD  (1,2,  3,4)  , nR3 

1 HI  a DT 

H2  a nl*0.500 
H3  a Hi  *2 .3  00 
H6  a HI  /6.  300 
g(M3  ) a 0.00 
A a 0.5  00 
GO  TO  5 

2 A a 0.2928922138134525 
GC  TO  8 

3 A a 1.7071067311365475 
GC  TO  5 

4 X('«3)  a X( '<3  I *H6«X00T  ( M3  )-0('^3  ) /3.C0 
RKL3E3  a 1. 

IF  ( ICONT.SC.NEO)  RKLDE2a2. 

3:  TO  6 

5 X(M3)  3 X(M3)  ,.A*(  0TaxDOT(M3)-3(  M2)  ) 

Q(M3J  > H3*  AaxOOT  (M3)  *(  1 .00-3  .00*A  )*Q(M3) 

RKL3E2  a 1, 

6 RETUON 
END 


FUNCTI  ON  CCPLX  (P,Z,  G,CR  IVE,X2,0MG,  NR4) 
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implicit  realms  (A-H,  0-Z) 

REALMS  THA0UT(3  001»3)  t XC  ATA ( 3 OQi  , 3 ) , C OCT  ( i ) 

REALM'S  CMG(  1 ) » P ( 1 I , Z ( 1 ) ,G  ( 1 ) , D R1  VE  ( 1 ) ,X  2 ( 2 ST  , X2D  CT  ( 25  ) 
DIMENSION  IOU.T(3) 

COMMON  THAOUT.XOA  TA»T  ,OT  ,TF  , 

CNOUT,  I IN,  MV  ,M3,  ICGNTtNEQ,  I SKIP,  ITF,  lOUT 
QMG0  0T(M3)  = X2  (M3) 

SS  = RKLD£3(0MG,0.mgDQT,NRA) 

X2D0T(M3)  = -2.^P(M3)*2(M3)*X2(M2)-Z(M3)»»2*0MG(M3)> 
1G(M3)*0RIVE( M3  ) 

SSS  = RKLDE4{X2,X2D0T,NRA) 

CCPLX  = SSS 

RETURN 

END 


1 

2 

3 

W 

5 

6 


FUNCTION  RKLDE4  (X,XDCT,NRA) 

IMPLICIT  REALMS  (A-H,0-Z) 

REAL*8  thAOUT  (3  001,3  ) , X OAT  A ( 300  1 , 2 ) 

REAL*8  X(l) , XOOT(l),  CC(25) 

OIMENS  ION  IOUT( 3) 

COMMON  THAOU'^tX  OATA,T,OT,TF, 

CNOUT, 1 1 N,NV  ,M3,IC0NT,NEa,I SKI P ,ITF, lOUT 
GO  TO  ( 1,2,  3,4)  , NR4 
HI  = OT 
H2  » Hl*0,500 
H3  = HI  *2.000 
H6  = Hi/6.000 
0C<M3)  3 0.00 
A = 0 .500 
GC  TO  5 

A = 0.2828922138134525 
GC  TO  5 

A = 1.7071067811865475 
GO  TO  5 

X(M3)  = X(  m3  )+h6*XD0T(M3)-QC(M3)/2.00 
RKL0E4  = 1. 

IF  ( ICONT.EO.NSg)  RKLDE4=2. 

GC  TO  6 

X(m3)  = X(  M3)+A*  ( OT*XDOT(M3  )-QC  ( m3  ) ) 

gC(M3)  =»  H3*A*X00T(  M3  )♦(  1.00-3.0  0*A)*QC(  M3) 
RKL054  = 1 . 

RET'JRM 

END 


FUNCTION  KE(C) 

EVALUATES  IMOLICIT  CONSTRAINTS  FOR  8GX  PLX 

REAL*3  C(25) 

KEaO 
R ETURN 
END 


FUNCT ION  FE (C ) 

THE  P UNCTION  MINIMIZED  SV  BOX  PLX 
implicit  R£AL*8  ( A-H,0-2  ) 

REAL *3  THAr)UT(300l,3)  ,XDATa(300:  ,3),C(25),W(3) 
CGMMnNj  THA0'JT,X0ATA,T  ,OT,'^F, 
*N0UT,IIN,NV,M3,ICaNT,NcC,ISKlP,  IT  = , I'^UT 
COMMON  /REGl/  W,IPM 
DC  10  I=1,NCUT 
THAOUTIl  ,I  ) =0.00 

10  continue 

OIFF*O.CO 

PI=0.00 

CALL  PLAN’CC) 

T =07 


A 
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DG  11  1=1, ITF 
00  12  J=l,N0UT 

DIFF=XD4TA(  I , J )- THA  GUT  ( I , J ) 
GO  TOd  ,2,3,^  ),  IPM 

1 PI  =PI*-W  (J)*(01FF*2) 

GO  TO  12 

2 PlaPI+W  (J  »♦  ( QIFF*=*‘2 
GO  TO  12 

3 P I = PI+W  (J)*0ABS(0  IFF  > 

GO  TO  12 

A PI  »PI  + W(  J)*DA3S(OIFF  )*T 

12  continue 

T»T+OT 
11  CONTINUE 
FE=P  I 
RETURN 
END 


SUBROUTINE  PPLT (NPPLT , ICP,  1 1) 

AUTOSCALES  AND  CALLS  FOR  PRINTR  PLOTS 

INPLI Cl T REALMS  (A-H,C-ZI 
R£AL*A.  XX(900),rY(900)  ,WW(  900) 

REAL«<i-  TX(A),TY(A-) 

REAL  *3  THA0UT(300  1,3)  , XDATA  ( 3 001  ,3  ) 
OIMENSICN  IQUT(3) 

COMMON  THA0UT,X3ATA,T  ,OT,TF, 

CNOUT,  I I N ,NV  ,M3  tlCONT,  NEC,I  SK I P , I TF  , I OUT 
OC  1 1=1,900 
XX(  I ) = 0. 

YY( I ) = 0. 

1 WW( I ) = 0. 

T=0.  DO 

TSTEP  = 5.00«0T 

J = 0 

a iGx=  0.00 

3IGY  = 0.D0 

SMLX=O.DO 

SMLY=0.00 

OC  2 I=1,I0P,5 

J = J + 1 

XXf  J ) * T 

YY(  J ) =XCATA  ( I,  I I ) 

WW(  J)  =THAOUT(I  ,11) 

X = T 

XO  = YY  (J) 

TH  = WW( J) 

Y«AX*0MAX1 (XD,TH) 

YMIN  = 0MIM  (X0,TH) 

XvaX  = 0MAX1 ( BIGX  ,X  ) 

I F (B  IGX.lt. X)  BIGX=X 
IF  OIGY.LT.YMAX)  3IGY  = Y«AX 
IF  (SMLY.GT.YMIN)  SMLY*YMIN 
T » t+tSTEP 

2 CONTINUE 
TX( 

T X( 

TX(3)  » SMLX 

TX(^)  » BIGX 
TY(n  * BIGY 
TY(5)  » SMLY 

TYO)  * 0. 

TY(V)  = 0. 

WRITE  (6,3)  BIGX,  BIGY  , SMLY 
CALL  olOTo  (Tx,TY,-i„l  I 
CALL  PLCTP  (XX,  WW,NPPLT  ,2) 

CALL  PLOTP  ( XX,  YY  ,iNPPLT,3) 

MBIT=  (6,A) 

0(  I5p.  £0.4500)  WRIT5(6,5) 
pET’JRN 


0. 

0. 
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3 format  ( 2X,  • BIGX=»  SElS.TfEXt’SIGYs  ',615. 7, 2X, 
l'S«LY=  ’,615.71 

4 CQR'^AK  //,2  X,'  SYST6M  R6SPaNR6  FOR  FRCBL6f' ') 

5 FOR'-^AT  ( //,2X,.'ST0P  AT  900  GRAPH  POINTS’) 

6  NO 


SUBR0UTIN6  PIC  (NOOLT , ICP,  1 1) 

AUTO  SCALES  ANO  CALLS  FOR  VERSA  TEC  PLOTS 

IMPLICIT  REAL*3  (A-H,C-Z) 

PEAL^’V  XX(  900  ),  YY  (900)  ,V,W(  900) 

RtAL-^TXIA)  ,Ty(  4) 

REAL*8  THAO  UT  (3001,3)  , XDATA  ( 3 001 , 2 ) , TI  TLE  ( 1 2 ) 

REAL  *3LABC/’  • / 

real  *4LABA/'  a '/,LABC/'  0 •/ 

OIMENSION  I0UT(3) 

COMMON  THAO'JT,XOATA,T  ,OT,TF,  ' 

*NOUT,I  IN,NV  ,M3,IC0NT,  NEC  , I SKI  = , IT  F » I OUT 
READ  (5,2)  TITLE 
T=0.00 

TSTEP  =5.00*aT 
J = 0 

3IGX  = 0.00 

3IGY  = o.no 

S'^LX  = 0.00 

S.^LY  = O.DO 

00  1 Ial,I0P,5 

J = J>1 

XX  (J)  a T 

YY(J)  = XOATA(I,in 

WW(J  )=’'hAOUT(  I,  II  ) 

X =«  T 

XO  * YY( J) 

T h s WW ( J ) 

YMAX  = CMAXl(XO,Th) 

YMIN  a OMINl ( XO,TH) 

XMAX  = CMAXK  0IGX,X  ) 

IF  (BIGX.LT.X)  3IGX=X 
IF  (3IGY.LT.YMAX)  3IGY  = YMAX 
IF  (SmlY.GT.YMIN ) SMLY=YMIN 
T = T>TSTE= 

1 continue 

’X( 1 ) = 0. 

TX(2)  = 0. 

"X(3)  =»  SMLX 
TX(4)  = 31GX 
TY(1)  = 3IGY 
TY(2)  = SMLY 
TY(3 ) » 0. 

TY(4)  3 C. 

CALL  OOAW  (4,TX,’'Y,  1,  1,LASC,TITLE  , C,G,  C,  C,  C,  0,3 ,3,0,L 
CALL  DRAW  (NPPLT,XX,WW  ,2,0,LASA,TrLS,0,J,0,0,a,0,8,a 
»0  ,L  ) 

CALL  DRAW  (NPOLT,  XX,  YY,3,0,L  A30,T  rtE,  C,  C,  0,0,0, 0, 3, 3 
*0,L) 

IF  (L.NE.O)  WRITE  (5,3)  L 
RETURN 

2 =ORMAT  ( 5A3  ) 

3 FORMA'-  (//,'  GRAPH  NOT  COMPLETED.  OUTPUT  CODE  = ',12 
END 


SUBROUTINE  30XPLX  ( MV , N AV , NPR  , N’’ Z , R Z , XS  , I P , 3U , 3 L , YMN  , 
•I  cR) 


SUBROUTINE  BOXPLX  (CATEGORY  HO) 

PURPOSE 


A 
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30XPLX  IS  A SUBROUTINE  USED  TQ  SOLVE  THE  dogblem 
OF  LOCATING  A MINIMUM  (OR  MAXIMUM)  OF  AN  ARBITRARY 
OBJECTIVE  FUNCTION  SUBJECT  TO  ARBITRARY  EXPLICIT 
AND/OR  implicit  CONSTRAINTS  BY  THE  COMPLEX  METHOD 
OF  M.j.  acx.  EXPLICIT  CONSTRAINTS  .ARE  DEFINED  AS 
UPPER  AND  LOWER  BOUNDS  ON  THE  INDEPENDENT  VARIABL- 
ES. implicit  constraints  M4Y  35  AQBITRARY  FUNCTIONS 
OF  THE  VARIABLES.  TWO  FUNCTION  SUBPROGRAMS  TC 
EVALUATE  THE  OBJECTIVE  FUNCTION  AND  IMPLICIT  CON- 
STRAINTSt  RESPECTI VELY,  MUST  RE  SUPPLIED  BY  THE 
USER  (SEE  example  BELOW).  3GXPLX  ALSO  HAS  THE  OPT- 
ION TO  PERFORM  INTEGER  PROGRAMMING,  WHERE  THE  VAL- 
UES OF  THE  INDEPENDENT  VARIABLES  ARE  RESTRICTED  TC 
INTEGERS  . 

USAGE 

CALL  BCXPLX  ( NV , N AV , NPR ,NT A , R , XS,  IP, XU, XL , YMN , lER  I 
OESCPIPTION  OF  PARAMETERS 

NV  AN  INTEGER  INPUT  DEFINING  THE  NUMBER  OF  INDE- 
PENDENT VARIABLES  OF  THE  OBJECTIVE  FUNCTION 
TO  BE  MINIMIZED.  NOTE;  m^xImUM  NV  + NAV  IS 
PPESENTLY  50.  MAXIMUM  NV  IS  25.  IF  THESE 
LIMITS  MUST  BE  EXCEEDED,  PUNCH  A SOURCE  DECK 
IN  the  usual  manner,  AND  CHANGE  THE  CIMEM- 

s ION  statements. 

NAV  AN  INTEGER  INPUT  DEFINING  "^HE  NUMBER  OF  AUXI- 
LIARY variables  the  user  WISHES  TO  DEFINE  FOR 
HIS  OWN  convenience.  '^YPICALLY  HE  M4Y  WISH 
TO  DEFINE  THE  VALUE  OF  EACH  IMPLICIT  CON- 
STRAINT FUNCTION  AS  AN  AUXILIARY  VARIABLE.  IF 
THIS  IS  DONE,  THE  CP"^!  ONAL  OUTPUT  FEATURE 
OF  BOXPLX  CAN  BE  USED  TO  OBSERVE  THE  VALUES 
OF  those  constraints  AS  THE  SOLUTION  PRO- 
GRESSES. AUXILIARY  VARIABLES,  IF  USED,  SHOULD 
BE  EVALJATED  IN  FUNCTION  KE  (DEFINED  BELOW). 
NAV  MAY  Sc  iERG. 

NPF  INPUT  INTEGER  CONTROLLING  THE  FREQUENCY  OF 
■UV'touT  DESIRED  FCR  DIAGNOSTIC  PURPOSES.  1 = 

.NOR  .LE.  J,  NO  OUTOU’  WILL  BE  PRODUCED  BY 
BOXPLX.  OTncRWlSe,  THE  CURRENT  CQmolEX  DF 
K*P*NV  VERTICES  AND  THEIR  CcN-RGIO  WILL  BE 
■•U^PUT  AFTER  EACH  '(PR  PERMISSIBLE  "^RiALS.  THE 
NUMBER  OF  total  TRIALS,  NUMBER  OF  FEASIBLE 
trials,  number  of  function  EVALUATIONS  ANC 
NUMBER  DF  implicit  GGN  STRAIN'  EVA  LJA’’ IONS  ARE 
I ncluoed  I n the  output  . 

AODI'^  IONALLY,  ( rtHEN  NPR  ,GT.  0)  THE  SAVE  IN- 
FORMATION WILL  BE  OUTPUT: 

1)  IR  THE  INITIAL  POINT  IS  NOT  FEASIBLE. 

2)  AFTER  The  FIRST  COMPLETE  COMPLEX  IS 
GENER  ATEO  . 

3)  IF  A FEASIBLE  VERTEX  CANNQT  3E  FOLNO  AT 
SOME  trial. 

A.)  IF  THE  OBJECTIVE  VALUE  CF  A VE.R'EX  CAN- 
NOT BE  -IADS  NO-LON  GER-wPRST  . 

5)  1=  THc  LIMIT  ON  tri^LS  (NTA)  IS  REACHED 
AMO , 

6)  when  the  objective  fijncTION  HAS  BEEN  UN- 
CHA.NGEO  FOR  2*NV  tpulS,  INDICATING  A 
LOCAL  MINIMUM  HAS  BEEN  FOUND. 

1=  THE  USER  WISHES  TO  'RACE  "SHE  PRCGRESS  OF 
A SOLUTION,  A CHOI SE  OF  NPC=  25,  5C  OR  100  IS 
RECCMMENOED. 

NTA  INTEGER  INPUT  CF  LI  MI'  ON  -^HE  NUMBER  DF 
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C C 
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R 

XS 

IP 

XU 

XL 

YMN 

lER 


EXAMOL 


TRIALS  ALLOWED  IN  ^hE  CALC  UL  A^^ION  . IF  THE 
USER  INPUTS  NTA  .L=.  Ot  A CSFAULT  VALUE  OF 
2000  IS  USED.  WHEN  THIS  LINIT  IS  REACHED 
CONTROL  RETURNS  TO  THE  CALLING  PRGGRAH  WITH 
THE  BEST  attained  OBJECTIVE  FUNCTION  VALUE  IN 
Y'^N,  AND  THE  BEST  ATTAINED  SOLUTION  OCINT  IN 
XS. 

A REAL  number  INPUT  TO  DEFINE  THE  FIRST  RAN- 
DOM number  USED  IN  DEVELOPING  THE  INI- 
TIAL COMPLEX  OF  Z’i'NV  VERTICES. 

(0.  GT.  R .LT.  1.)  IF  R IS  NOT  WITHIN  THESE 
BOUNDS,  IT  WILL  BE  REPLACED  BY  1./3.  . 

INPUT  REAL  ARRAY  DIMENSIONED  AT  LEAST  NV+NAV. 
THE  FIRST  NV  MUST  CONTAIN  A FEASIBLE  ORIGIN 
FOR  STARTING  THE  CALCULA""!  GN.  THE  LAST  NAV 
NEED  NOT  BE  INITIALIZED.  UPON  RETURN  FROM 
BOXPLX,  the  FIRST  NV  ELEMENTS  OF  THE  ARRAY 
CONTAIN  THE  COORDINATES  OF  ‘I’HE  MINIMUM  OB- 
JECTIVE FUNCTION,  AND  THE  REMAINING  NAV 
(NAV  .GE.  0»  contain  THE  VALUES  OF  T|-e 
CORRESPONDING  AUXILIARY  VARIABLES. 

INTEGER  INPUT  FOR  OPTIONAL  IN'^EGER  RRO- 
GRAMMING.  if  IP=1,  the  values  of  THE  INDE- 
PENDENT VARIABLES  WILL  BE  REPLACED  WITH 
integer  VALUES  (STILL  STORED  AS  R£AL*4). 

A REAL  ARRAY  DIMENSIONED  AT  LEAST  NV  INPUT- 
TING THE  UPPER  BOUND  ON  EACH  INOE=£NDENT 
VARIABLE,  (EACH  EXPLICI'  CON  STR  A I N^^  ) . INPUT 
VALUES  ARE  SLIGHTLY  ALTERED  BY  BOXPLX. 

A PEAL  ARRAY  DIMENSIONED  AT  LEAST  NV  INPUT- 
TING THE  LOWER  BOUND  ON  EACH  INOEOENDENT 
VARIABLE,  (EACH  EXPLICIT  CCNSTRAINT).  NOTE: 
FOR  BOTH  XU  AND  XL  CHOOSE  REASONABLE  VALUES 
IF  NONE  ARE  GIVEN,  NOT  VALLES  WHICH  ARE  MAG- 
NITUDES ABOVE  OR  BELOW  THE  EXPECTED  SOLUTION. 
INPUT  values  are  SLIGHTLY  ALTERED  EY  BOXPLX. 

THIS  OUTPUT  IS  THE  VALUE  (REAL*4)  OF  the  OB- 
JECTIVE FUNCTION,  CORRESPONDING  TO  THE  SOLU- 
TION POINT  OUTPUT  IN  XS. 

INTEGER  ERROR  RETURN.  tq  pc  INTERROGATED 
UPON  RETURN  FROM  BOXPLX.  lER  WILL  BE  ONE  OF 
the  FOLLOWING: 

=-l  CANNOT  FIND  FEASIBLE  VERTEX  CR  FEAS- 
IBLE CENTRHID  AT  THE  START  OR  A RE- 
START (SEE  'mcthOC  BELOW). 

=0  FUNCTION  VALUE  UNCHANGED  FOR  'N' 

trials.  (WHERE  N = 6»NV-t-ia)  THIS  IS 
The  NORMAL  RETURN  PARAMETER. 

= 1 CANNOT  DE'^ELOP  FEASIBLE  VERTEX. 

=2  CANNGT  DEVELOP  A NO-LCNGEP-wOOST  VER- 
TEX. 

= 3 LIMIT  TRIALS  REACHED.  (NTA  EX- 
CEEDED) 

NQTC:  valid  RESULTS  MAY  BE  RETURNED  IN 

ANY  OF  THE  ABOVE  CASES. 

OF  USAGE 


C 

r 

c 

c 

c 

c 

r 

C 


c 

c 

c 
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W 

c 

c 

c 

c 

r 

C 

c 

c 
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C 

C 
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r 

C 

C 
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c 

r 

C 

c 

c 

c 

r 

'r 

c 

c 

c 

c 

c 


THIS  example  minimizes  the  QBjeC'IVE  FUNCTION  C 

SHOWN  IN  THE  external  FUNCTION  F£(X).  THERE  ARE  C 
TWO  independent  variables  X (1  ) C X(2I,  AND  TWO  I M-  C 
PLICIT  CONSTRAINT  FUNCTIONS  X(3)  & X(4)  WHICH  ARE  C 
EVALUATED  AS  AUXILIARY  VARIABLES  (SEE  EXTERNAL  C 

= UNCTION  KE(X)  ).  C 


J 
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0I»^EMSI0N  XS  (4)  » XU(  2)  f XL(  21 


STARTING  GUESS 
XS(l)  = 1.0 
XS{2)  = 0.5 
UPPER  LIMITS 
XU(1)  = 6.0 
XU(2  ) = 6.  0 
LOWER  LIMITS 
XL(l  ) = 0.0 
XL(2)  a 0.0 


R = 9./15. 
NTA  = 5000 
NOR  = 50 
NAV  = 2 
NV  = 2 
10  = 0 


CALL  80XPLX  ( N V , NA  V , N PR  ,NTA  , R » XS  . I P,  XU  t XL  . YMN  . I ER  ) 
WRIT  6 (6,1)  ((XS(  I ),  1=1,4)  ,YMN,  lER  ) 

1 FORMAT! ////  ,'  THE  PCINT  IS  LOCATED  AT  (XS<I»=) 
*4(E13.7,5X)  ,//,  • A.NO  THE  FUNCTION  VALUE  IS 
♦ E13.7  , • lER  = ',15) 

STOP 

END 


PUNCTICN  KS  (X) 

EVALUATE  CONSTRAINTS,  SET  KE=0  IF  NC  IMPLICIT 
CONSTRAINT  IS  VIOLATED,  CR  SET  i<E=l  IF  ANY  IMPLICIT 
CONSTRAINT  IS  VIOLATED. 

dimension  X(4) 

XI  = X ( 1 ) 

X2  = X(2) 

X E = 0 

X (3  ) = XI  1.732051»X2 

IF  (X(3)  .LT.  0.  .OR.  X(3)  . GT . 6.)  GO  TO  1 
X(4)  = Xl/1. 732051  -X2 
IF  (X ( 4)  .G£.  0 . ) RETURN 

1 KE  = 1 
RETURN 
END 


FUNCTION  FE(X) 

DIMENSION  X(4) 

THIS  IS  THE  09JECTIVE  ^UNCTION. 

FE=  -(X(2)*»3  * (9.-(  X(  l)-3.  )**2  ) / (46.T3533)  ) 

R ET'JP  N 
E NO 

METHOD 

THE  C CmolEX  METHOD  IS  AM  EXTENSION  AND  ADAPTI CN  OF 
'HE  simplex  method  C."  LINEAR  ^RPORAMMING.  START- 
ING WITH  ANY  ONE  FEASISLE  POINT  IN  N-CIMENSIDN 
SPACE  A ’’COMPLEX'*  OF  2*N  VERTICES  IS  CONStR'jC'ED 
3y  SELEC'^’ING  random  ROINTS  within  the  feasible 
REGION.  PCR  THIS  PURPOSE  N CCCPOINATES  ARE  FIRST 
randomly  chosen  WITHIN  THE  SPACE  BOUNDED  3Y  EX- 
PLICIT CONSTRAINTS.  THIS  DEFINES  A TRIAL  INITIAL 
VERTEX.  IT  IS  then  CHECKED  FOR  POSSIBLE  VIOLA'ION 
CF  implicit  CONSTRAINTS.  IF  ''NE  CR  mq^e  aRE  VIO- 
LATED, THE  trial  initial  VEO'^EX  IS  DISPLACED  nALF 
OF  ITS  distance  from  The  CENTRQIO  OF  PREVIOUSLY 
SELECTED  INITIAL  VERTICES.  IF  NECESSARY  THIS  CIS- 


69 


oooooof  ’oooo*  ‘To  >of  ■><1000,  It  >000000000<-i00<-,o0o00o00o0t“>o00o0o000o0o00o000000 


r«  >oo<''»<'*o<'*f  »i  '*oooo<">oor>0‘~»<^o<  icfooc  *riooo<  'oooo<''ooooo<  >f->i-»c'>oooooooooooof''oot'>o<'>oor)0 


placement  process  is  repeated  until  the  vertex  has 

BECOME  FEASIBLE.  IF  THIS  FAILS  TC  HAPPEN  AF^ER 
5*N-«-10  OISPLACEMENTSt  THE  SOLUTION  IS  ABANDONED. 
A^-^ER  EACH  VERTEX  IS  ACDEO  TO  THE  COMPLEX.  THE 
CURRENT  CENTROID  IS  CHECKED  FOR  FEASIBILITY.  IF 
IT  IS  INFEASIBLE,  THE  LAST  TRIAL  VERTEX  IS  ABAND- 
ONED AND  AN  effort  TC  GENERATE  AN  ALTERNATIVE 
TRIAL  VERTEX  IS  MADE.  IF  5^N+1C  VERTICES  ARE 
ABA iN DONE 0 CONSECUTIVELY,  THE  SOLUTION  IS  TER- 
MINATED. 

IF  AN  INITIAL  COMPLEX  IS  ESTABLISHED,  THE  BASIC 
COMPUTATION  LOOP  IS  INITIATED.  THESE  INSTRUCTIONS 
FIND  THE  CURRENT  WORST  VERTEX,  THAT  IS,  THE  VERTEX 
WITH  THE  LARGEST  CORRESPONDING  VALUE  FOR  THE  OB- 
JECTIVE FUNCTION,  AND  REPLACE  THAT  VERTEX  BY  ITS 
OVER-REFLECTION  THROUGH  THE  CENTROID  OF  ALL  OTHER 
VERTICES.  (IF  THE  VERTEX  TO  BE  REPLACED  IS  CON- 
SIDERED AS  A VECTOR  IN  N-SPACE,  ITS  OVER- 
REFLECTION IS  OPPOSITE  IN  DIRECTION,  INCREASED  IN 
length  by  the  factor  1.3,  AND  CCLLINEAR  WITH  THE 
REPLACED  VERTEX  AND  CENTROID  OF  ALL  OTHER 
VERTICES.  ) 

WHEN  AN  OVER-REFLECTION  IS  NOT  FEASIBLE  OR  REMAINS 
WORST,  IT  IS  CONSIOEREC  NQt-p ER M I S3  I 2L E ANO  IS 
DISOLACEO  halfway  TOWARD  the  CEN-"R0ID.  AFTER  FCUR 
SUCH  ATT£mpt5  are  MACE  UNSUCCESSFULLY,  EVERY  FIFTH 
ATTEMPT  IS  MADE  BY  REFLECTING  THE  OFFENDING  VERTEX 
THROUGH  the  PRESENT  BEST  VERTEX,  INSTEAD  OF 
THROUGH  the  CENTROID.  IF  5*N+10  DISPLACEMENTS  AND 
OVER-REFLECTIONS  OCCUR  WITHOUT  A SUCCESSFUL 
(PERMISSIBLE)  RESULT,  THE  CURR£^'T  3EST  VERTEX  IS 
TAKEN  AS  AN  INITIAL  FEASIBLE  POINT  FOR  A RESTART 
RUN  OF  THE  COMPLETE  PROCESS.  RESTARTING  IS  ALSO 
U'NOERTAKEN  when  5*NV+10  CONSECUTIVE  TRIALS  HAVE 
BEEN  MADE  WITH  NO  SIGM=ICAN'^  CHANGE  IN  the  VALUE 
OF  THE  OBJECTIVE  FUNCTION.  IN  ALL  CASES,  RESTART- 
ING IS  INHIBITED  IF  THE  LAST  RESTART  rj  n ^JQr  pqq- 
DUCE  A SIGNIFICANT  IMPROVEMENT  IN  th£  VINIMUm 
ATTAINED. 

IT  IS  RECOMMENDED  THAT  THE  USER  READ  The  REFER- 
ENCE FOR  FURTHER  USEFUL  INFORMATION.  IT  SHOULD  BE 
NOTED  that  THE  ALGORITHM  DEFINED  THERE  HAS  SEEN 
ALTERED  TO  FIND  THE  CONSTRAINED  RATHER 

THAN  THE  MAXIMUM. 


REMARKS 


THE  INTEGER  PROGRAMMING  OPTION  WAS  ADDED  TQ  THIS 
PROGRAM  AS  SUGGESTED  IN  R5FEO£NCE  (2).  A MIXED 
INTEGER/CONTINUOUS  VARIABLE  VERSION  CF  BCXPLX 
WOULD  PE  EASY  TO  CREATE  BY  DECLARING  ”1?''  TO  BE  AN 
ARRAY  OF  NV  CONTROL  VARIABLES  WHERE  IP(I)=1  WOULD 
INDICATE  that  the  I-TH  VARIABLE  IS  TO  BE  CONFINED 
TO  INTEGER  values.  EACH  STATEMENT  0^  THE  FORM 
•1=  (IP  .EC.  1)'  ETC.  WOULO  then  NEED  tq  3?  AL- 
TERED TO  'IF  (IP(I)  .£3.  1)'  ETC.,  WHERE  "XE  SUB- 
SCRIPT IS  APPROPRIATELY  CHOSEN.  NORMALLY,  XU  AND 
XL  VALUES  ARE  ALTERED  TO  BE  AN  EPSILON  ‘WITHIN* 
actual  values  declared  BY  THE  USER.  THIS  ADJUST- 
*^ENT  IS  NOT  made  when  IP*1. 


NOTE:  NO  NON-LINEAR  PROGRAMMING  ALGORITHM  CAN 

GUARANTEE  ThAT  th£  ANSWER  =OUMO  IS  THE  GLOBAL 
MINIMUM,  rather  Than  just  a local  'minimum.  how- 
ever, ACCORDING  TO  R£F.  2,  THE  COMPLEX  METHOD  HAS 
AN  advantage  in  THAT  IT  TENDS  "^D  =IND  THE  GLOBAL 
MINIMUM  MOPE  FREDUENTLY  THAN  MANY  OTHER  NON- 
LINEAR programming  ALGORITH'^S. 
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IT  SHCULD  BE  NOTED  Th4T  ThE  AUXILIARY  VARIABLE 
FEATURE  CAN  ALSO  BE  USED  DEAL  WITH  PRG8LEMS 
CONTAINING  EQUALITY  CONSTRAINTS.  ANY  EOUAL^IY 
CONSTRAINT  IMPLIES  THAT  A GIVEN  VARIABLE  IS  NOT 
TRULY  INDEPENDENT.  THEREFORE,  IN  GENERAL,  ONE 
VARIABLE  INVOLVED  IN  AN  EQUALITY  CONSTRAINT  CAN 
BE  RENUMBERED  FROM  THE  SET  OF  NV  INDEPENDENT 
VARIABLES  AND  ADDED  TO  THE  SET  OF  NAV  AUXILIARY 
VARIABLES.  THIS  USUALLY  INVOLVES  RENUMBERING 
THE  INDEPENDENT  VARIABLES  OF  the  GIVEN  PROBLEM. 

SUBROUTINES  AND  FUNCTIONS  REQUIREC 

SUBROUTINE  ’BOUT*  ANC  FUNCTION  'FBV  ARE  INTEGRAL 
OARTS  OF  THE  BOXPLX  PACKAGE. 

TWO  FUNCTIONS  MUST  EE  SUPPLIED  BY  THE  USER.  THE 
FIRST,  KE(X),  IS  USED  TO  EVALUATE  ""HE  I-PLICIT 
CONSTRAINTS.  SET  K E=0  AT  THE  BEGINNING  OF  THE 

function,  then  evaluate  the  I'^olICIT  constraints. 

IN  THE  example  ABOVE,  THE  FIRST  CONSTRAINT,  X O)  , 
MUST  BE  WITHIN  THE  RANGE  (0.  .L  E . X(3)  .LE.  6). 
THE  SECONC  CONSTRAINT  X(-V),  MUST  BE  .GS.  0.  . IF 
EITHER  CONSTRAINT  IS  NOT  WITHIN  THESE  3CUN0S, 
CONTROL  IS  TRANSFERRED  TO  STATEMENT  1,  ANC  KE  IS 
SET  TO  '’I"  AND  CONTROL  IS  RETURNED  TO  BOXPLX. 

THE  SECOND  FUNCTION  ThE  USER  must  oROVIDE  EVALU- 
ATES THE  OBJECTIVE  FUNCTION.  IT  IS  CALLED  FE(X) 
AS  SHOWN  IN  THE  EXAMPLE  ABOVE,  AND  FE  MUST  35  SET 
TO  THE  VALUE  OF  THE  OBJECTIVE  pUNCT’ION  CCRRES- 
POMDING  TO  CURRENT  VALUES  OF  THE  NV  INDEPENDENT 
VARIABLES  IN  ARRAY  »X  • . 
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implicit  REALMS  (A-H,0-Z» 

REAL* 3 V(50  ,50J  ,FUN( 30)  ,SUM(25 ) ,CEN  (25J  ,XS (NV  ) , 3U{ NV  ) , 
*8L(NV) 

KV  » 5 
E®  * l.D-6 
NTA  = 2000 

1=  (NTZ.GT.O)  NTA  = NTZ 

■5  = 37 

IFtR.LE.O.DO .OR .R .G£ .1  .CO )R  = 1 .DO/ 2 .00 
N vT  3 N V+NA  V 

total  VARS,  EXPLICIT  PLUS  I mol  Id  T 

NT  * 0 

CURRENT  TRIAL  NO. 

NFT  3 0 

CURRENT  NO.  OF  PERMISSIBLE  tpialS 
NTFS  3 0 
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CURRENT  NG.  GF  TIMES  F HAS  BEEN  ALMOST  UNCHANGED 


CHECK  FEASIBILITY  OF  START  POINT 

DC  I = ltNV 
VT  = XS(I ) 

IF  I BL  ( I)  .LE.VT  ) GO  TO  1 
II  » -I 
VT  a 3H  n 
GO  TO  2 

1 IF  (BUdl.GE.VT  ) GO  TC  3 
I I a I 

VT  a 0U(  n 

2 IF  (NPR.GT.O)  WRITE  (6,49)  II 

3 Vdtl)  a VT 
C EN  ( I ) =»  VT 

I F ( IP. EQ.1 ) go  to  4 

BL(I ) = BL (1) +0MAX1( EP, EP*0ASS(BL  (I  ) ) ) 

3U( I )=8U( I )-OMAXl ( EP,  EP  + CABS ( BU ( 1)  ) ) 

4 SLM(  n a VT 


NCE  a 1 

number  of  constraint  evaluations 

I * 1 

IF  (K5 ( V(l,l) ). Sa.O)  GC  TO  5 
I F (NOR  ,LE  .0)  GO  TO  12 
WRITE  (6,50) 

GO  TO  12 

5 NFS  a 1 

NUMBER  OF  VERTICES  ( K)  a 2 '^IMES  NG.  CF  VARIABLES. 

K a 2 ♦NV 

number  of  displacements  allowed. 

NLIM  = S^NV+IO 

NUMBER  0*=  CONSECUTIVE  TRIALS  wl^^H  U^CHANGEC  Ft  TO 
TERMINATE. 

NCT  a NLIM+NV 
A LOHA  =1  . 3D  0 
FK  a K 
FKM=FK-1  .00 
BETAaALPHA^-  l.DO 

INSURE  SEED  OF  RANDOM  NUMBER  GENE'^ATOR  IS  OOC. 

i Qq  - R *1 , 3 7 

IF  (MO  D(  IQR  ,2  ) .S3  .J  ) I3RaI3o,.101 

SET  UP  initial  VER'^ICES 
FUN  (I  ) a PE(V( 1, 1 ) ) 

Ymn  .=  FUNd  ) 

6 FI  a l.oc 

FUNOLD  a F'Ji\(l) 

DO  15  Ia2,K 
FI  a FIwl.DO 
LIMT  a 0 

7 LIMT  a lIMT,-! 

END  CALCULATION  IF  FEASIBLE  CcNTRQIC  CANNC*^  BE  FOUND. 
IF  ( L Imt  .GE  .NLI  d GO  TO  11 

DC  3 J=1  ,NV 

C RANDOM  NUMBER  GENERATOR  (RAN DU) 

I CR  a I Gpa65539 

IF  (IQR.LT.O)  I3R  a IOR+21 47433647  + 1 
RCX  a IQR 

RCX  a R3X*.  465661  30-9 

V(J,I)  a 3L(J)+RQXM(au(J)-3L(  J?  ) 

iF({'/(  J,I)+.5  00  ) .3T.2147433647  . )WRITE(  6,  100) 
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IF  (ID.EQ.l  )V(J,n  = IDINT(V(  Jt  I ) + .5D0) 

3 CONTINUE 

OC  10  L=1»NHM 
NC5  = NCE+1 

IF  (KE(V(1,  I)  ).EQ.O)  GO  TO  13 
00  9 J=l,NV 

VT  = (V(J,I)KEN(  J)  »*.5D0 

IF(  ( VT  + .5D0  J.GT  .2 14748  3647  . JW  R IT  P ( 6, 100  ) 

IF  (IP.EO.l)  VT^IOINT  (VT+.500) 

V ( J, n = VT 
9 CONTINUE 

10  CONTINUE 

11  IF  (NPR  .LE. 0»  GO  TO  12 
WRITE  (6»51)  I 

CALL  BOUT  (NT,NPT  ,NFE,NCE,NV,NVT,V,I,PJN,CEN,  I) 

12  lER  = -1 


13  00  14  J = 1,MV 

SUM(J)  = SUM(J) +V(J,  I) 

14  C£N(  J)  a SUM(  J)  /FI 

TRY  TO  ASSURE  FEASIBLE  CENTROID  FOR  STARTING. 
NCS  a NCE+1 

IF  (KE(CEN).EQ.0>  GO  TO  60 
suM(j)  = suM(j»  -v( J,  n 
GO  TO  T 
60  NFE  = MFE+1 

FUN(  I ) a FE  (V(l  , I n 

15  CONTINUE 


END  CF  LOOP  SETTING  OF  INITIAL  COMPLEX. 

I P ( NPR.LE.  0)  GO  TO  17 

CALL  BOUT  (NT  ,NPT  ,NF£,  NCE,  NV,N  VT,  V,K,FUN»CEN,  0) 

PIND  THE  WORST  VERTEX,  THE  'J'TH. 

J = 1 


00  16  1=2, K 

IF  (FUN  (J  ) .GE.FiJ.N(  I ) ) 30  16 

J = I 

16  CCNTIMUE 

BASIC  LOOP,  eliminate  EACH  WORST  VERTEX  IN  TURN.  IT 
MUST  BECOME  NO  LONGER  WORST, NOT  MERELY  IMPROVED.  FIND 
NEX  T-t-i-wORS’  VERTEX, the  ’JN'^h  :)N|g. 

17  JN  a 1 

IF  (J.EQ.l)  JN  a 2 


OC  13  I al  ,K 

IF  ( I .EQ.J  ) GO  TO  13 

IF  (FJN  (JN  ) .GE.FUNI  n ) GO  tq  la 

JN  a I 

13  CONTINUE 

LIMTaNUMSER  OF  MOVES  DURING  THIS  TRIAL  TOW  ARC  TrE 
CENTROID  DUE  TO  FUNCTION  VALUE. 

LT'T  a 1 

COMPUTE  CENTRQIO  AND  OVER  REFLECT  WCRS^  VERTEX. 


00  19  I *1  ,NV 
VT  a V(  I ,J) 

SUM(  I ) a 3UM(  I J-VT 

CEN(  I ) a SUM(i»  /PKM 

VT  * 5ETA*CEN(  I )-ALPHA*VT 

1 F(  (VT  + .5  00  ).GT  .2147  433  647  . )W  RITE  (6,100) 
IF  (IP.EO.l)  VT  arDINT( VT+.5D0  ) 


I 
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C 

C INSURE  ^HE  EXOLICIT  CONSTRAINTS  ARE  0RSERV6D. 

IR  VUfJ)  = OMAXK  DMINHVT»EU  II ) ) , BL  U ) ) 

C 

NT  = NT+1 

CHECK  PGR  IHOLICIT  CONSTRAINT  VIOLA-^ICN. 

20  OC  25  N*1»NLIM 
NCE  = NCE+1 

IF  «E{V(1,J)).E3  .0)  GO  TG  26 
EVERY  'KV'TH  TI  Me  »0  VER-REFCECT  THE  OFFENDING  VERTEX 

through  the  best  vertex. 

IF  (MOOINtKV)  .NE.OI  GO  TC  22 
CALL  FBV  (K, FUN, Ml 

DC  21  I*1,NV 

VT  a BETA*V(I,M)-ALPHA*V(I,J)- 
IF(  (VT  ,..500  ) .GT.2147^a2647.  IWR  ITE(  6,  100) 

IF  (IP.EQ.l)  VT  =IDI  NT  (VT^-.5D0  ) 

21  V(I,J)  * 0MAX1(QMIN1(VT,8U(  I)  ) ,0L(  1)  ) 

GC  TG  24 

CONSTRAINT  VIOLATION:  MOVE  NEW  POINT  TOWARD  CENTROID. 

22  OG  23  I*1,NV 
VT  a (CEN(  I ) >V  ( I,  J I )*.5  00 
I = ( { VT-t-,  500)  .GT.214748  364T.  ) WR  IT=  (6 ,100  ) 

IF  (lo.EQ.l)  VT  = lOIN T( VT+.500) 

V (I  , J ) * 7T 

23  CONTINUE 

NT  = MT+1 

25  CONTINUE 

I ER  = 1 

CANNOT  GET  FEASIBLE  VERTEX  BY  MOVING  TOWARD  CEN'RCID, 

OR  SY  OVER-REFLECTING  THRU  THE  BEST  VERTEX. 

1=  (NPR.L5.0)  GO  TO  42 
'WRITE  ( 6,  52)  NT,  J 

CALL  BOUT  (NT,NPT ,NFE,MC£,NV,NVT, v,<,FUN,CcN, J) 

GC  TO  42. 

FEASIBLE  vertex  FOUND,  EVALUATE  THE  OBJECTIVE  ‘=UNCTION. 

26  NFE  = NFE'*’! 

F'JMTRY  « FE(  V(  1,  J ) ) 

TEST  TO  SEE  IF  FUNCTION  VALUE  HAS  NCt  CHANGED. 

AFO  =0ABS  ( PJNTR  Y-FUNOLD  ) 

AWXaOMAXl  I DABS  (cP*FUNCLC  ) , £3  ) 

ACTIVATE  THE  FOLLOWING  TWO  STATEMENTS  FOR  DIAGNOSTIC 
PURPOSES  ONLY. 

WRITE  (6, 90)  J,  AFO,  AM  X,  FUNTR  Y,FUMOLC,  PiJN  < J ) ,FLN(  JN)  , 

♦ NTFS  ,N 

99  FCRMAT  (IX, 13,6=15.7,215) 

IF  ( AFO.  GT  .AMx)  GO  TO  27 
NTFS  = NTFS 

IF  (NTFS.LT.NCT)  GO  TO  28 
IcR  = 0 

IF  (NPR.Lr.O)  GO  TO  42 
WRITE  (6,53)  K 

CALL  BOUT  (Nt,npT ,NFE,NCE,MV,NVT, V,K,FUN,C£N, C) 

GO  TO  42 

27  NTFS  » 0 

IS  THE  N£w  VER*^EX  NO  LCNG£R  WORST? 

28  IF  (FUNTRY  .LT.FUN  (JN)  ) GO  TO  34 
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C ^RIAL  VERTEX  IS  STILL  WORST;  ADJUST  TOWARD  CENTROID. 

C EVERY  'KV’TH  TIME,  OVER-RE^LEC^  Tl-6  CE^ENOING  VER'^EX 

C THROUGH  the  BEST  VERTEX. 

LIMT  = LIMT+1 

IF  (MOO (LIMT,KV) .NE.O ) GO  tq  30 
CALL  F3V  (K,FUN,M) 

C 

00  29  I =1,NV 

VT  » BETA*V ( I,M l-AL?HA*V( I ,J) 

IF(  (VH-.5C0  ).GT. 2 147483  647.  )W  RITE  (6»100) 

IF  (IP.EQ.l)  VT  al  OINT(  VT+.500  ) 

29  V(I,J)  s DMAXK  DM  INK  VT,BU(  I ) » ,BL(  I ) ) 

C 

GO  TO  32 
C 

30  00  31  I»1,NV 

VT  = (CEN(  I)+V(  I,  J ) )<.500 

IF(  (VT+.300J  .GT. 2147483647.  )WR  ITE(6,100  » 

1=  (IP.EQ.l)  VT  alDI  NT(  VT+.500) 

V ( I,  J ) =»  VT 

21  ccrriNUc 

r 

32  IF  (LIMT.LT.NLIM)  GO  TO  33 

: CANNOT  MAKE  THE  ' J*  TH  VERTEX  NO  LONGER  WORST  3Y 

C DISPLACING  TOWARD  T rt£  CENTROID  OR  BY  OVER-REFLECTING 
C THROUGH  the  BEST  VER'^EX. 
lER  = 2 

IF  (NPR  .L£.  0)  GO  TO  42 
WRIT?  (6,5  2)  NT,  J 

CALL  BOUT  ( N‘!’,NPT,NFE  ,NCE,N  V,N  VT  , V ,K,FUN  ,C£N,  J) 

GO  TO  42 
23  NT  3 NT'fl 
GO  TO  20 

SUCCESS;  W£  HAVE  A REPLACEMENT  FOR  VERTEX  J. 

34  FUN(J)  » FUNTRY 
FUNOLD  = FUNTRY 
NPT  = NPT+1 

EVERY  IQO'Th  permissible  TR  I AL  , R ECO  voiJTc  CENTROID 
summation  to  avoid  CREEPING  ERROR. 

IF  (Mno(  vipT,  IJO  ) .NE  .0  ) GO  TO  37 

DC  36  I=1,NV 
SUM(  I ) = 0.00 

DO  35  N »1  ,K 

25  5UM(  I ) 3 $UM(  I)<-V(I,N) 

.:EN(  I ) * SUNK  I ) /fk 

36  continue 

1 C = 0 
GO  TO  39 

37  DO  38  I 3l  ,NV 
28  S;JM{I  ) 3 SUM(  I)+V(  I,J) 

LC  = J 

30  IF  (NPR  .LE  .0  ) GO  TO  40 

IF  (MOD(  NPT  ,NPR  ) . NE.O  ) GC  t^  <vo 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,C5N,  LC  ) 

HAS  THE  MAXIMUM  NUMBER  CF  TRIALS  BEEN  REACHED  WIThCU"^ 
convergence?  if  not,  30  TO  NEW  T?IAL. 

40  IF  (NT.GE.NTA)  go  to  41 

NEXT-TO-WORST  VERTEX  NOW  BECOMES  WORST. 

J » JN 
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GC  TO  17 
lER  = 3 

I F (NOR  ,GT  .0  ) WR  ITE  (6,54) 


COLLECTOR  OQI^^T■FaR  ALL  ENDINGS. 

1)  CANNOT  develop  FEASIBLE  VERTEX. 

2)  CANNOT  DEVELOP  A NC-LCMGEP-WOOST  VERTEX 

3)  FUNCTION  value  UNCHANGED  FOR  K TRIALS. 

4)  LIMIT  ON  TRIALS  REACHED. 

5)  CANNOT  PINO  FEASIBLE  VERTEX  AT  START. 

42  CONTINUE 


I PR 

i'r 


lER 

lER 

lER 


2 

0 

3 

-I 


FIND  BEST  VERTEX. 

CALL  FBV  (K,FJN,M) 

IF  ( IER.GE.3  ) GO  TO  44 


restart  IF  THIS  SOLUTION  IS  SIGNIFICANTLY  BETTER  thaN 
THE  PREVIOUS, OR  IF  THIS  IS  THE  FIRST  TRY. 

IF  (NPR.LE.O)  go  to  43 
-JRITE  (6,55)  ( M ,Y  MN,FUN(  M)  ) 

43  IF  (FUN( M) .gE.YMN  ) GO  TC  47 

IF(DA3S(FUN( M)-yNN) .LE.DMAX: (EP, EF*YMN) ) GO  TO  47 


GIVE  I"!*  ANO'^HER  TRY  UNLESS  LIMIT  ON  TRIALS  REACHED. 

A4  Y'-N  = FUN(M) 

FUN(l)  » FJN(M) 

DO  45  I »1 ,NV 
C £N(  I ) = V ( I,  M) 

SUM (I ) = V ( I,M) 

45  v( 1 ,1 ) * v(  I ,m; 

DC  ^6  1=1,  NVT 

46  XSd  ) = V(  I ,M) 

IF  (IER.lt. 31  GO  TO  6 

47  1=  ( NOR.LE.  0)  GO  TO  48 

CALL  BOUT  (NT.NOT  ,iNFE,NCE,NV,NV“,Y,K,FUN,V(  1,M)  ,-l) 

WRITE  (6.56)  F^J^^(v,) 

48  RETURN 
C 

49  format  ( 'OINCEX  AND  DIRECTION  OF  '^U'^LYING  VARIABLE', 

* ' AT  START'  ,I  5) 

50  FORMA"' ( *0  Imol  ic  IT  CONSTRAINT  VIOLATED  AT  START.', 

*0  X,  • DEAD  END.  • ) 

51  FORMAT!  • OCANMQT  FIND  FE  A SI  BLE  ' , 14  , 

*'TH  VERTEX  '''R  CENTROIC  AT  START.'  ) 

52  =GRMAT  ( ' OA'''  TRIAL  ',14,'  CANNOT  FIND  FEASIBLE  VERTEX', 
*'WHICH  IS  NO  LONGER  WOR  ST ' , I 4 , 1 5 X , • RE  ST  AR  T FROM  BEST', 
«•  VERT  EX  ,'  ) 

53  =0RMAT  (AOHOFUNCTIQN  HAS  BEEN  AL^CST  UNCHANGED  FOR  15, 
*7H  TR I ALS) 

54  FCRMA"'  (2‘^HOLIMIT  CN  TRIALS  EXCEECEC.  ) 

55  FORMAT!  • 38  = ST  VERTEX  IS  NO . ' , I 3 , ' CLD  min  WAS 
*615.7,  ' NEW  min  is  ' ,E13.7) 

56  =GRMAT  I'OMIN  OBJECTIVE  “JNCtIQN  IS  ',£15.7) 

100  FORMAT ( 'O'  , 'TRUNCATI ON  ERROR  IN  BGXPLX') 

END 


SUBROUTINE  FBV  (i<,pUN,m) 
IMPLICIT  REAL»3  (A-M,0-Z) 

= 6AL*8  FUN(  50) 

Mai 

00  I 1=2, K 

IF  (FIJN(  M)  .lE  .FUN  ( I ) ) GO  TO  1 
M = I 

1 COMT INU  E 
3£"'URN 
END 
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5UBPGUTIME  BOUT  ( N T, N P T ,NF E ,NC E ,N V ,NVT  , V , K tFN  , C , I K ) 
implicit  peaces  (4-H,0-ZJ 
REALMS  V(  50,50)  ,Frv|(50)  ,C(25) 

XRITE  (6,4)  MT,iMPT,NF£,NC£ 

DC  1 I =1,K 

WRITE  (6,5)  FN(  I)  ,(  V(  J,I  ) ,J*1  ,NV) 

IF  (NVT.LE.NV)  GO  TO  1 
MVP  * NV+1 

WRITE  (6,6)  ( V(  J,  I ) , J=NVP,NVT) 

1 CCNTINUE 
C 

IF  ( IK  .ME.O  ) GO  TO  2 
C 

WRITE  (6,7)  (C(I  ) ,I  =l  ,NV) 

RETURN 

2 IF  ( IK.GE.O  ) GO  TO  3 
WRITE  ( 6,3)  (C(  I)  ,I  =1  ,NV) 

RETURN 

3 WRITE  (6,R)  IK,  (C(  I ) , I = 1,NV  ) 

R ETURN 

r 


4 FCRMATCONO.  TOTAL  TRIALS  = ',I5,4X, 

♦ ’NO.  FEASIBLE  TRIALS  = •,I3,4X, 

♦ •NO.  FUNCTION  EVALUATIONS  = ' , I 5,  AX, 

♦'NO.  constraint  EVALUATIONS  = ',15/, 

♦'  C function  V ALUE'  ,6X, 

♦ • INDEPENDENT  V ARI A8L £S / C £P ENDEN T QR  IMPLICIT’, 
♦'C0NSTRAIN‘'’S’  ) 

5 format  ( 1h  ,Eia.7,2X,7E14.7/(llX,7El-i.7)  ) 

6 FORMAT  (21X  , 7E1  <i-.7  ) 

7 =ORMAT  ( lOHCCEiNTROID  1 IX  ,7E14 . 7/ ( 21 X , 7E 14 . 7 ) ) 

8 FORMA-^  CO  REST  V£R  TE  X’ ,7X,  7E  14.  7/ ( 2iX  , TEJA.  7)  ) 

9 FORMAT  ('0C5NTR3ID  LESS  VX  • , 1 2 , 2X  , 7E1 4.  7/ ( f IX,  7 5 14 . 7 ) ) 
END 

C 

C A /*  CARD  GOES  HERE 

C 

//GO.SYSIN  00  ♦ 
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